!!$!!!!!!!!!!!!!!!!!!!!!!!! !!$!!!! Copyright 2002 Lorenzo Giada and Matteo Marsili !!$!!!! lgiada@mpikg-golm.mpg.de marsili@sissa.it !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!$ This program is free software; you can redistribute it and/or modify !!$ it under the terms of the GNU General Public License as published by !!$ the Free Software Foundation; either version 2 of the License, or !!$ (at your option) any later version. !!$ !!$ This program is distributed in the hope that it will be useful, !!$ but WITHOUT ANY WARRANTY; without even the implied warranty of !!$ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!$ GNU General Public License for more details. !!$ !!$ You should have received a copy of the GNU General Public License !!$ along with this program; if not, write to the Free Software !!$ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module mod_grpdet integer :: N, idum, Nc, mfl, clnm integer, allocatable, dimension(:) :: a,m,ags,m_gs integer, allocatable, dimension(:,:) :: indx real :: egs real, allocatable, dimension(:) :: x,e,x_gs,e_gs real, allocatable, dimension(:,:) :: C character(11) :: filfin,fildat,filout character(10) :: filgs character(9) :: filord character(3) :: filxxx character(7) :: filin character(18), allocatable, dimension(:) :: cells end module mod_grpdet !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! program group_fit !!$ calculates the gs deterministically: !!$ 1. merges two clusters with maximal likelihood; !!$ 2. deterministic single step moves; !!$ 3. go to 1. until there is only 1 cluster !!$ REQUIRES the files: !!$ grpdet.par !!$ covariance matrix xxx.cij !!$ names of series (optional) xxx.names (written as 2 columns: number, name) !!$ WRITES the files: !!$ xxx.ent.det the energy every sweep of the lattice, and other run-time info !!$ xxx.smile the energy after step 1. and after step 2. !!$ xxx.cnf.mrg minimal energy configuration from merging moves ONLY !!$ xxx.cnf.det the ground state configuration after all single moves use mod_grpdet implicit none integer :: Nk,ugh integer,allocatable,dimension(:) :: m0,a0,indi,a00 integer :: i,j,k,l,ik,nt,cell integer :: mi_new,mf_new,ai,af,nupd integer :: m0_new,m0_min,i0_min,j0_min real, allocatable, dimension(:) :: ee real :: xi_new,xf_new,ei_new,ef_new,de,etot real :: e0_min,e0_new,C0_new,C0_min,e00 real, allocatable, dimension(:) :: e0 real, allocatable, dimension(:,:) :: C0 !!$ READ data and name of files from: grpdet.dat call ini(ugh) if (ugh.eq.1) goto 998 allocate(e0(N),m0(N),a0(N),C0(N,N),ee(N),a00(N)) egs=10000.0 e00=1000.0 !!$ initial condition UNIFORM do i=1,N a(i)=i end do call cond if (mfl.ne.1) open(20,file=filout,status='unknown') a0=a m0=m e0=e C0=C Nc=N 50 e0_min=10000.0 !!$ select 2 clusters to be merged a=a0 call cond !THIS is important to reconstruct the array indx do i=1,N-1 if (m0(i).gt.0) then do j=i+1,N if (m0(j).gt.0) then C0_new=C0(i,i)+C0(i,j)+C0(j,j) m0_new=m0(i)+m0(j) call energy(m0_new,C0_new,e0_new) if (e0_new-e0(i)-e0(j).lt.e0_min) then C0_min=C0_new e0_min=e0_new-e0(i)-e0(j) i0_min=i j0_min=j end if end if end do end if end do !!$ MERGE C0(i0_min,i0_min)=C0_min do k=1,m0(j0_min) a0(indx(j0_min,k))=i0_min indx(i,m0(i0_min)+k)=indx(j0_min,k) end do m0(i0_min)=m0(i0_min)+m0(j0_min) call energy(m0(i0_min),C0_min,e0(i0_min)) !not necessarily in this way do k=1,N if (k.ne.i0_min.and.k.ne.j0_min) then C0(i0_min,k)=C0(i0_min,k)+C0(j0_min,k) C0(k,i0_min)=C0(i0_min,k) end if end do !!$ remove the index j0min m0(j0_min)=0 e0(j0_min)=0.0 C0(j0_min,j0_min)=0.0 do k=1,N C0(j0_min,k)=0.0 C0(k,j0_min)=0.0 end do if (mfl.eq.1) then Nc=Nc-1 open(90,file=filxxx//'.smile',position='append') write(90,*) Nc, sum(e0)/real(N) close(90) if (sum(e0).lt.e00) then a00=a0 e00=sum(e0) end if goto 777 ! skip single data moves end if !!$ Starting point a=a0 call cond call gs Nc=Nc-1 !!$ccc CICLE of single moves. call cic_temp Nk=0 ! number of clusters now do i=1,N if (m(i).gt.0) Nk=Nk+1 end do open(90,file=filxxx//'.smile',position='append') write(90,*) Nc, sum(e0)/real(N), Nk, sum(e)/real(N) close(90) if (sum(e0).lt.e00) then a00=a0 e00=sum(e0) end if 777 if (Nc.gt.1) go to 50 !!$ done: close file with energy if (mfl.ne.1) close(20) !!$ writes GS from merging open(99,file=fildat) do i=1,N write(99,*) a00(i) end do close(99) a=a00 call cond k=0 do i=1,N if(m(i).gt.0) k=k+1 end do write(6,*) 'minimal energy from cluster-merging',e00/float(N), 'in ', k,'clusters' if (mfl.eq.1) goto 998 !finished !!$ Saves the GS in file filfin (cnf_...) open(12,file=filfin) do i=1,N write(12,*) ags(i) end do close(12) k=0 do i=1,N if(m_gs(i).gt.0) k=k+1 end do write(6,*) 'there are ',k,'clusters, and the GS energy is ',sum(e_gs)/real(N) 998 continue end program group_fit !!$cccccccCCCCCCCCcccccccCCCCCCcccccccCCCCCccccccc subroutine sweep(nupd) !!$ SWEEP of the lattice use mod_grpdet implicit none integer :: nupd integer :: i,k,j,l,ai,af,ncls,ik,afmin integer, dimension(0:1) :: mi_new,mf_new real :: de,dxf,demin real, dimension(0:1) :: xi_new,ei_new,xf_new,ef_new do i=1,N l=i !sequential sweep ai=a(l) ik=0 !label for trick demin=1. do af=1,N !try all moves = all other clusters if (m(af).eq.0) go to 88 ! This is FORBIDDEN in KMINES ! if (m(af).eq.0.and.m(ai).eq.1) go to 88 ! simple relabelling: abort if (af.eq.ai) go to 88 ! no change of cluster: abort if (m(ai).eq.1) then mi_new(ik)=0 xi_new(ik)=0.0 ei_new(ik)=0.0 else xi_new(ik)=x(ai) do k=1,m(ai) j=indx(ai,k) xi_new(ik)=xi_new(ik)-C(l,j) end do mi_new(ik)=m(ai)-1 call energy(mi_new(ik),xi_new(ik),ei_new(ik)) end if mf_new(ik)=m(af)+1 dxf=0.0 do k=1,m(af) j=indx(af,k) dxf=dxf+C(l,j) end do xf_new(ik)=x(af)+1+dxf call energy(mf_new(ik),xf_new(ik),ef_new(ik)) de=ei_new(ik)+ef_new(ik)-e(ai)-e(af) if (de.lt.demin) then afmin=af demin=de ik=mod(ik+1,2) end if 88 end do ik=mod(ik+1,2) if (demin.lt.0) then ! accept move af=afmin !!$ update index: 1. remove the index i from cluster ai do k=1,m(ai)-1 if (indx(ai,k).eq.l) indx(ai,k)=indx(ai,m(ai)) end do indx(af,mf_new(ik))=l ! 2. add new index !!$ update cluster variables m(ai)=mi_new(ik) x(ai)=xi_new(ik) e(ai)=ei_new(ik) m(af)=mf_new(ik) x(af)=xf_new(ik) e(af)=ef_new(ik) !!$ update Potts variable a(l)=af nupd=nupd+1 end if end do !end lattice sweep return end subroutine sweep !!$cccccccCCCCCCCCccccccCCCCCccccccccCCCCcccCCCCCc subroutine cic_temp use mod_grpdet implicit none integer :: nupd integer :: i,k,j,l,ai,af,t,t0,ncls real :: mav !!$ccc CICLE IN TEMPERATURE t=0 10 nupd=0 call sweep(nupd) t=t+1 mav=0.0 ncls=0 do i=1,N if (m(i).gt.0) ncls=ncls+1 mav=mav+m(i) call energy(m(i),x(i),e(i)) ! this is just for fun end do !!$cccccc WRITES initial # of clusters, number of sweeps , avg. # updates, energy, !!$cccccc # clusters, avg. occupation of clusters in file filout (xxx.ent) write(20,999) Nc,t,float(nupd)/(float(N)),sum(e)/real(N),ncls,& &mav/ncls 999 format(2I5,f10.4,f12.6,i5,f12.6) if (nupd.gt.0) go to 10 call gs return end subroutine cic_temp !!$cccccccCCCCCCCCCcccccccccCCCCCCCccccccCCCCCCCccccc subroutine cond use mod_grpdet implicit none integer i,j,k !!$ Calculates parameters for a given cluster configuration, !!$ which is in the array a(N) m=0 x=0.0 e=0.0 do i=1,N m(a(i))=m(a(i))+1 end do do i=1,N k=1 do j=1,N if (a(j).eq.i) then indx(i,k)=j k=k+1 end if end do end do do i=1,N if (m(i).eq.1) x(i)=1. if (m(i).gt.1) then do j=1,m(i) do k=j,m(i) !! this is important, because C is already twice x(i)=x(i)+C(indx(i,j),indx(i,k)) end do end do end if end do do i=1,N if (m(i).gt.1) then call energy(m(i),x(i),e(i)) end if end do return end subroutine cond !!$ccccccCCCCCCccccccccccCCCCCCCcccccccccCCCCCCCCCccccc subroutine ini(ugh) use mod_grpdet implicit none integer :: i,j,ugh !!$cccc READ parameters open(11,file='grpdet.par',status='old',iostat=ugh) if(ugh.ne.0) then write(6,*) 'configuration file grpdet.par is missing. will STOP' goto 999 end if read(11,*) N read(11,*) filxxx read(11,*) clnm read(11,*) mfl close(11) allocate(C(N,N),indx(N,N)) allocate(a(N),m(N),ags(N),x(N),e(N),x_gs(N),m_gs(N),e_gs(N),cells(N)) filin(1:3)=filxxx filin(4:7)='.cij' filout(1:3)=filxxx filout(4:11)='.ent.det' filfin(1:3)=filxxx filfin(4:11)='.cnf.det' fildat(1:3)=filxxx fildat(4:11)='.cnf.mrg' filord(1:3)=filxxx filord(4:9)='.names' !!$ LOAD CORRELATION MATRIX !!$ OBSERVE: multiply by 2 is correct open(11,file=filin,status='old') do i=1,N-1 do j=i+1,N read(11,*) C(i,j) C(i,j)=2*C(i,j) C(j,i)=C(i,j) end do C(i,i)=1.0 end do C(N,N)=1.0 close(11) !!$ load name of series if(clnm.eq.1) then open (12,file=filord,status='old', iostat=ugh) if(ugh.ne.0) then write(6,*) 'cannot find file with name of series' write(6,*) 'please check variable clnm in input file grpdet.par' write(6,*) 'will continue without names' clnm=0 goto 888 end if do i=1,N read(12,*) j,cells(i) end do end if 888 continue ugh=0 return 999 ugh=1 return end subroutine ini !!$cccccccCCCCCCCccccccccCCCCCcccccCCCCCCCCC subroutine gs !!$ updates the ground-state use mod_grpdet implicit none integer :: k, ik real :: etot etot=sum(e) if (etot.lt.egs) then ags=a x_gs=x m_gs=m e_gs=e egs=etot end if return end subroutine gs !!$cccccccccCCCCCCCCCCCCCCcccccccccccccCCCCCCCc subroutine energy(m,x,e) !!$ calculate the energy of a cluster implicit none integer :: m real :: x,e,chi if (m.le.1) then e=0.0 !!$cCccc controllo che c_s