| 1 | c**************************************************************
|
|---|
| 2 | c
|
|---|
| 3 | c This file contains the subroutines: partem_p
|
|---|
| 4 | C Compared to the version in the main distribution, this
|
|---|
| 5 | C routine doesn't write the rmsd nor native contacts to the time
|
|---|
| 6 | C series.
|
|---|
| 7 | c
|
|---|
| 8 | c Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 9 | c Shura Hayryan, Chin-Ku Hu
|
|---|
| 10 | c Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 11 | c Jan H. Meinke, Sandipan Mohanty
|
|---|
| 12 | c
|
|---|
| 13 | c **************************************************************
|
|---|
| 14 |
|
|---|
| 15 | subroutine partem_p(num_rep, nequi, nswp, nmes, nsave, newsta,
|
|---|
| 16 | & switch, rep_id, partem_comm)
|
|---|
| 17 | C
|
|---|
| 18 | C PURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM
|
|---|
| 19 | C ON PARALLEL COMPUTERS USING MPI
|
|---|
| 20 | C
|
|---|
| 21 | C switch: Choses the starting configuration:
|
|---|
| 22 | C -1 - stretched configuration
|
|---|
| 23 | C 0 - don't change anything
|
|---|
| 24 | C 1 - random start configuration
|
|---|
| 25 | C
|
|---|
| 26 | c CALLS: addang,contacts,energy,hbond,helix,iendst,metropolis,
|
|---|
| 27 | c outvar,(rand),rgyr
|
|---|
| 28 | C
|
|---|
| 29 | include '../INCL.H'
|
|---|
| 30 | include '../INCP.H'
|
|---|
| 31 | include '../incl_lund.h'
|
|---|
| 32 | include 'mpif.h'
|
|---|
| 33 |
|
|---|
| 34 | logical newsta
|
|---|
| 35 | integer switch, partem_comm, rep_id, nsave
|
|---|
| 36 |
|
|---|
| 37 | external can_weight
|
|---|
| 38 |
|
|---|
| 39 | C nequi: number of Monte Carlo sweeps for thermalization
|
|---|
| 40 | C nswp: number of Monte Carlo sweeps
|
|---|
| 41 | C nmes: number of Monte Carlo sweeps between measurments
|
|---|
| 42 | C newsta: .true. for new simulations, .false. for re-start
|
|---|
| 43 |
|
|---|
| 44 | dimension eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC),
|
|---|
| 45 | & inode(MAX_PROC), geavm(MAX_PROC), gsph(MAX_PROC)
|
|---|
| 46 | double precision pbe(MAX_PROC),yol(MAX_PROC),acy(MAX_PROC),
|
|---|
| 47 | & acy1(MAX_PROC),acx1(MAX_PROC),
|
|---|
| 48 | & rgyrp(MAX_PROC), eol0,rgyp,acz0
|
|---|
| 49 |
|
|---|
| 50 | double precision e_min, e_minp(MAX_PROC), e_minpt(MAX_PROC)
|
|---|
| 51 | integer h_max, h_maxp(MAX_PROC)
|
|---|
| 52 | c Order of replica exchange
|
|---|
| 53 | integer odd
|
|---|
| 54 | ! Counter to keep random number generators in sync
|
|---|
| 55 | integer randomCount
|
|---|
| 56 |
|
|---|
| 57 | c Collect partial energies. Only the root writes to disk. We have to
|
|---|
| 58 | c collect the information from the different replicas and provide
|
|---|
| 59 | c arrays to store them.
|
|---|
| 60 | c eyslr storage array for solvent energy
|
|---|
| 61 | c eyelp - " - coulomb energy
|
|---|
| 62 | c eyvwp - " - van-der-Waals energy
|
|---|
| 63 | c eyhbp - " - hydrogen bonding energy
|
|---|
| 64 | c eysmi - " - intermolecular interaction energy
|
|---|
| 65 | double precision eyslr(MAX_PROC)
|
|---|
| 66 | double precision eyelp(MAX_PROC),eyvwp(MAX_PROC),eyhbp(MAX_PROC),
|
|---|
| 67 | & eyvrp(MAX_PROC),eysmip(MAX_PROC)
|
|---|
| 68 | c Collect information about accessible surface and van-der-Waals volume
|
|---|
| 69 | c asap storage array for solvent accessible surface
|
|---|
| 70 | c vdvolp storage array for van-der-Waals volume
|
|---|
| 71 | double precision asap(MAX_PROC), vdvolp(MAX_PROC)
|
|---|
| 72 |
|
|---|
| 73 | integer nhelp(MAX_PROC),nbetp(MAX_PROC), mhbp(MAX_PROC),
|
|---|
| 74 | & nctotp(MAX_PROC)
|
|---|
| 75 | integer imhbp(MAX_PROC)
|
|---|
| 76 | character*80 filebase, fileNameMP, tbase0,tbase1
|
|---|
| 77 | c frame frame number for writing configurations
|
|---|
| 78 | c trackID configuration that should be tracked and written out
|
|---|
| 79 | c dir direction in random walk
|
|---|
| 80 | c -1 - visited highest temperature last
|
|---|
| 81 | c 1 - visited lowest temperature last
|
|---|
| 82 | c 0 - haven't visited the boundaries yet.
|
|---|
| 83 | c dirp storage array for directions.
|
|---|
| 84 | integer frame, trackID, dir
|
|---|
| 85 | integer dirp(MAX_PROC)
|
|---|
| 86 |
|
|---|
| 87 | frame = ifrrm
|
|---|
| 88 | trackID = 1
|
|---|
| 89 | odd = 1
|
|---|
| 90 | write (*,*) 'Starting parallel tempering.'
|
|---|
| 91 | write (*,*) 'parameters, ',switch,newsta,nmes,nswp,nmes,
|
|---|
| 92 | & rep_id, num_rep, partem_comm, myrank
|
|---|
| 93 | call flush(6)
|
|---|
| 94 | C
|
|---|
| 95 | c
|
|---|
| 96 | C File with temperatures
|
|---|
| 97 | open(11,file='temperatures_abeta',status='old')
|
|---|
| 98 |
|
|---|
| 99 | tbase0='trj_00000'
|
|---|
| 100 | open(18,file=fileNameMP(tbase0,5,9,rep_id),status='unknown')
|
|---|
| 101 | if (rep_id.eq.0.and.myrank.eq.0) then
|
|---|
| 102 | c File with time series of simulation
|
|---|
| 103 | open(14,file='ts.d',status='unknown')
|
|---|
| 104 | endif
|
|---|
| 105 |
|
|---|
| 106 | C READ IN TEMPERATURES
|
|---|
| 107 | do i=1,num_rep
|
|---|
| 108 | read(11,*) j,temp
|
|---|
| 109 | pbe(j) = 1.0d0/( temp * 1.98773d-3 )
|
|---|
| 110 | end do
|
|---|
| 111 | close(11)
|
|---|
| 112 |
|
|---|
| 113 | c nresi: number of residues
|
|---|
| 114 | nresi=irsml2(1)-irsml1(1)+1
|
|---|
| 115 | C
|
|---|
| 116 | C Initialize variables
|
|---|
| 117 | do i=1,num_rep
|
|---|
| 118 | acx1(i) = 0.0d0
|
|---|
| 119 | acy(i) = 0.0d0
|
|---|
| 120 | eavm(i) = 0.0d0
|
|---|
| 121 | sph(i) = 0.0d0
|
|---|
| 122 | geavm(i) =0.0d0
|
|---|
| 123 | gsph(i) = 0.0d0
|
|---|
| 124 | e_minp(i) = 1.0d15
|
|---|
| 125 | h_maxp(i) = 0
|
|---|
| 126 | dirp(i) = 0
|
|---|
| 127 | end do
|
|---|
| 128 | dirp(1) = 1
|
|---|
| 129 | dirp(num_rep) = -1
|
|---|
| 130 | e_min = 1.0d15
|
|---|
| 131 | h_max = 0
|
|---|
| 132 | dir = dirp(rep_id + 1)
|
|---|
| 133 |
|
|---|
| 134 | c _________________________________ Initialize Variables
|
|---|
| 135 | if(newsta) then
|
|---|
| 136 | iold=0
|
|---|
| 137 | do i=1,num_rep
|
|---|
| 138 | inode(i) = i
|
|---|
| 139 | intem(i) = i
|
|---|
| 140 | end do
|
|---|
| 141 | c _________________________________ initialize starting configuration
|
|---|
| 142 | if (switch.ne.0) then
|
|---|
| 143 | do i=1,nvr
|
|---|
| 144 | iv=idvr(i) ! provides index of non-fixed variable
|
|---|
| 145 | if (switch.gt.0) then
|
|---|
| 146 | dv=axvr(i)*(grnd()-0.5)
|
|---|
| 147 | vr=addang(pi,dv)
|
|---|
| 148 | else
|
|---|
| 149 | vr = pi
|
|---|
| 150 | endif
|
|---|
| 151 | vlvr(iv)=vr
|
|---|
| 152 | enddo
|
|---|
| 153 | endif
|
|---|
| 154 | else
|
|---|
| 155 | if(rep_id.eq.0.and.myrank.eq.0) then
|
|---|
| 156 | open(13,file='par_R.in', status='unknown')
|
|---|
| 157 | read(13,*) iold
|
|---|
| 158 | do i=1,num_rep
|
|---|
| 159 | read(13,*) j,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
|
|---|
| 160 | end do
|
|---|
| 161 | jold=(iold/nmes)*num_rep
|
|---|
| 162 | rewind 14
|
|---|
| 163 | do i=1,jold
|
|---|
| 164 | read(14,*) idum1,idum2,idum3,dummy
|
|---|
| 165 | & ,dummy, dummy, dummy
|
|---|
| 166 | & ,dummy, idum1, idum2, idum3
|
|---|
| 167 | call flush(6)
|
|---|
| 168 | end do
|
|---|
| 169 | close(13)
|
|---|
| 170 | end if
|
|---|
| 171 | CALL MPI_BCAST(IOLD,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
|
|---|
| 172 | CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
|
|---|
| 173 | CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
|
|---|
| 174 | CALL MPI_BCAST(YOL,num_rep,MPI_DOUBLE_PRECISION,0,
|
|---|
| 175 | # MPI_COMM_WORLD,IERR)
|
|---|
| 176 | CALL MPI_BCAST(E_MINP, num_rep, MPI_DOUBLE_PRECISION, 0,
|
|---|
| 177 | # MPI_COMM_WORLD, IERR)
|
|---|
| 178 | CALL MPI_BCAST(h_maxp,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
|
|---|
| 179 | $ IERR)
|
|---|
| 180 | end if
|
|---|
| 181 |
|
|---|
| 182 | BETA = pbe(inode(rep_id+1))
|
|---|
| 183 | e_min = e_minp(inode(rep_id+1))
|
|---|
| 184 | h_max = h_maxp(inode(rep_id+1))
|
|---|
| 185 | write (*,*) "E_min=",e_min," for ", intem(rep_id + 1)
|
|---|
| 186 | eol=energy()
|
|---|
| 187 | if(.not.newsta.and.abs(yol(rep_id + 1) - eol).gt.0.1) then
|
|---|
| 188 | write(*,*) rep_id, ' Warning: yol(rep_id).ne.eol:'
|
|---|
| 189 | write(*,*) rep_id, yol(rep_id + 1), eol
|
|---|
| 190 | endif
|
|---|
| 191 | C Start of simulation
|
|---|
| 192 | write (*,*) '[',rep_id, myrank, beta, partem_comm,
|
|---|
| 193 | & '] Energy before equilibration:', eol
|
|---|
| 194 | c =====================Equilibration by canonical Metropolis
|
|---|
| 195 | do nsw=1,nequi
|
|---|
| 196 | call metropolis(eol,acz,can_weight)
|
|---|
| 197 | end do
|
|---|
| 198 | CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
|
|---|
| 199 | write (*,*) '[',rep_id,'] Energy after equilibration:', eol
|
|---|
| 200 | call flush(6)
|
|---|
| 201 | C
|
|---|
| 202 | C======================Multiple Markov Chains
|
|---|
| 203 | acz = 0
|
|---|
| 204 | do nsw=1,nswp
|
|---|
| 205 | c------------First ordinary Metropolis
|
|---|
| 206 | call metropolis(eol,acz,can_weight)
|
|---|
| 207 | iold = iold + 1
|
|---|
| 208 | eol0 = eol
|
|---|
| 209 | if(mod(iold,nmes).eq.0) then
|
|---|
| 210 | if ((rep_id + 1).eq.trackID.and.myrank.eq.0) then
|
|---|
| 211 | frame = iold /nmes
|
|---|
| 212 | filebase = "frame_00000.pdb"
|
|---|
| 213 | call outpdb(0, fileNameMP(filebase, 7, 11, frame))
|
|---|
| 214 | endif
|
|---|
| 215 | acz0 = acz
|
|---|
| 216 | C Measure global radius of gyration
|
|---|
| 217 | call rgyr(0,rgy,ee)
|
|---|
| 218 | rgyp = rgy
|
|---|
| 219 | C Measure Helicity and Sheetness
|
|---|
| 220 | call helix(nhel,mhel,nbet,mbet)
|
|---|
| 221 | C Measure Number of hydrogen bonds
|
|---|
| 222 | mhb = 0
|
|---|
| 223 | do i = 1, ntlml
|
|---|
| 224 | call hbond(i,tmhb,-1)
|
|---|
| 225 | mhb = mhb + 1
|
|---|
| 226 | enddo
|
|---|
| 227 | call interhbond(imhb)
|
|---|
| 228 | C Measure total number of contacts (NCTOT) and number of
|
|---|
| 229 | C native contacts (NCNAT)
|
|---|
| 230 | call contacts(nctot,ncnat,dham)
|
|---|
| 231 | c Add tracking of lowest energy configuration
|
|---|
| 232 | if (eol.lt.e_min) then
|
|---|
| 233 | c Write out configuration
|
|---|
| 234 | i=rep_id+1
|
|---|
| 235 | j=inode(i)
|
|---|
| 236 | e_min = eol
|
|---|
| 237 | filebase = "c_emin_0000.pdb"
|
|---|
| 238 | call outpdb(0, fileNameMP(filebase, 8, 11, j))
|
|---|
| 239 | filebase = "c_emin_0000.var"
|
|---|
| 240 | call outvar(0, fileNameMP(filebase, 8, 11, j))
|
|---|
| 241 | filebase = "c_emin_0000.dat"
|
|---|
| 242 | open(15, file=fileNameMP(filebase, 8, 11, j),
|
|---|
| 243 | & status="unknown")
|
|---|
| 244 | ! write(15,'(i8,2i4,f6.2,2f8.2,5i8)') iold,i,j,pbe(i),
|
|---|
| 245 | write(15,*) iold,i,j,pbe(i),
|
|---|
| 246 | & eol, eysl, eyel, eyvw, eyhb, eyvr, eysmi,asa,
|
|---|
| 247 | & vdvol, rgy, nhel, nbet, mhb, imhb, nctot
|
|---|
| 248 | close(15)
|
|---|
| 249 | endif
|
|---|
| 250 | c Add tracking of configuration with larges hydrogen contents.
|
|---|
| 251 | if ((mhb + imhb).gt.h_max) then
|
|---|
| 252 | c Write out configuration
|
|---|
| 253 | i = rep_id + 1
|
|---|
| 254 | j = inode(i)
|
|---|
| 255 | h_max = mhb + imhb
|
|---|
| 256 | filebase = "c_hmax_0000.pdb"
|
|---|
| 257 | call outpdb(0,fileNameMP(filebase,8,11,j))
|
|---|
| 258 | filebase = "c_hmax_0000.var"
|
|---|
| 259 | call outvar(0,fileNameMP(filebase,8,11,j))
|
|---|
| 260 | filebase = "c_hmax_0000.dat"
|
|---|
| 261 | open(15, file=fileNameMP(filebase, 8, 11, j),
|
|---|
| 262 | & status="unknown")
|
|---|
| 263 | ! write(15,'(i8,2i4,f6.2,2f8.2,5i8)') iold,i,j,pbe(i),
|
|---|
| 264 | write(15,*) iold,i,j,pbe(i),
|
|---|
| 265 | & eol, eysl, eyel, eyvw, eyhb, eyvr, eysmi,asa,
|
|---|
| 266 | & vdvol, rgy, nhel, nbet, mhb, imhb, nctot
|
|---|
| 267 | close(15)
|
|---|
| 268 | endif
|
|---|
| 269 |
|
|---|
| 270 | C
|
|---|
| 271 | C--------------------Gather measurement data
|
|---|
| 272 | ! I only use the master node of each replica for data collection. The
|
|---|
| 273 | ! variable partem_comm provides the appropriate communicator.
|
|---|
| 274 | if (partem_comm.ne.MPI_COMM_NULL) then
|
|---|
| 275 | CALL MPI_GATHER(RGYP,1,MPI_DOUBLE_PRECISION,RGYRP,1,
|
|---|
| 276 | & MPI_DOUBLE_PRECISION, 0,partem_comm,IERR)
|
|---|
| 277 | CALL MPI_GATHER(NHEL,1,MPI_INTEGER,NHELP,1,MPI_INTEGER,
|
|---|
| 278 | & 0,partem_comm,IERR)
|
|---|
| 279 | CALL MPI_GATHER(NBET,1,MPI_INTEGER,NBETP,1,MPI_INTEGER,
|
|---|
| 280 | & 0,partem_comm,IERR)
|
|---|
| 281 | CALL MPI_GATHER(MHB,1,MPI_INTEGER,MHBP,1,MPI_INTEGER,
|
|---|
| 282 | & 0,partem_comm,IERR)
|
|---|
| 283 | CALL MPI_GATHER(iMHB,1,MPI_INTEGER,iMHBP,1,MPI_INTEGER,
|
|---|
| 284 | & 0,partem_comm,IERR)
|
|---|
| 285 | CALL MPI_GATHER(NCTOT,1,MPI_INTEGER,NCTOTP,1,MPI_INTEGER,
|
|---|
| 286 | & 0,partem_comm,IERR)
|
|---|
| 287 | CALL MPI_GATHER(dir,1,MPI_INTEGER,dirp,1,MPI_INTEGER,
|
|---|
| 288 | & 0,partem_comm,IERR)
|
|---|
| 289 | CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1,
|
|---|
| 290 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 291 | CALL MPI_GATHER(e_min,1,MPI_DOUBLE_PRECISION,e_minp,1,
|
|---|
| 292 | & MPI_DOUBLE_PRECISION,0, partem_comm,IERR)
|
|---|
| 293 | CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,
|
|---|
| 294 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 295 | CALL MPI_GATHER(eysl,1,MPI_DOUBLE_PRECISION,eyslr,1,
|
|---|
| 296 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 297 | CALL MPI_GATHER(eyel,1,MPI_DOUBLE_PRECISION,eyelp,1,
|
|---|
| 298 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 299 | CALL MPI_GATHER(eyvw,1,MPI_DOUBLE_PRECISION,eyvwp,1,
|
|---|
| 300 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 301 | CALL MPI_GATHER(eyhb,1,MPI_DOUBLE_PRECISION,eyhbp,1,
|
|---|
| 302 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 303 | CALL MPI_GATHER(eyvr,1,MPI_DOUBLE_PRECISION,eyvrp,1,
|
|---|
| 304 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 305 | CALL MPI_GATHER(eysmi,1,MPI_DOUBLE_PRECISION,eysmip,1,
|
|---|
| 306 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 307 | CALL MPI_GATHER(asa,1,MPI_DOUBLE_PRECISION,asap,1,
|
|---|
| 308 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 309 | CALL MPI_GATHER(vdvol,1,MPI_DOUBLE_PRECISION,vdvolp,1,
|
|---|
| 310 | & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 311 |
|
|---|
| 312 | c Write trajectory
|
|---|
| 313 | write (18,*) '@@@',iold,inode(rep_id+1)
|
|---|
| 314 | call outvbs(0,18)
|
|---|
| 315 | write (18,*) '###'
|
|---|
| 316 | ! call flush(18)
|
|---|
| 317 | c Write current configuration
|
|---|
| 318 | if ((mod(iold, nsave).eq.0)) then
|
|---|
| 319 | filebase = "conf_0000.var"
|
|---|
| 320 | call outvar(0, fileNameMP(filebase, 6, 9, rep_id+1))
|
|---|
| 321 | endif
|
|---|
| 322 | endif
|
|---|
| 323 |
|
|---|
| 324 | if(rep_id.eq.0.and.myrank.eq.0) then
|
|---|
| 325 | randomCount = 0
|
|---|
| 326 | c Update acceptance, temperature wise average of E and E^2 used to calculate
|
|---|
| 327 | c specific heat.
|
|---|
| 328 | do i=1,num_rep
|
|---|
| 329 | j=intem(i)
|
|---|
| 330 | acy(i)=0.0
|
|---|
| 331 | c Above: contents of acy1 are added to acy(i) a few lines down.
|
|---|
| 332 | c acy1(intem(i)) contains information received from the node at temperature
|
|---|
| 333 | c i, on how many updates have been accepted in node intem(i). Since acz
|
|---|
| 334 | c is not reset to 0 every cycle, acy(i) must be set to 0 here. Else, there
|
|---|
| 335 | c will be serious double counting and the values of acceptance printed
|
|---|
| 336 | c will be simply wrong.
|
|---|
| 337 | e_minpt(i)=e_minp(intem(i))
|
|---|
| 338 | end do
|
|---|
| 339 | do i=1, num_rep
|
|---|
| 340 | e_minp(i) = e_minpt(i)
|
|---|
| 341 | j=intem(i)
|
|---|
| 342 | acy(i)=acy(i)+acy1(j)
|
|---|
| 343 | eavm(i)= eavm(i)+yol(j)
|
|---|
| 344 | sph(i) = sph(i)+yol(j)*yol(j)
|
|---|
| 345 | enddo
|
|---|
| 346 |
|
|---|
| 347 |
|
|---|
| 348 | C Write measurements to the time series file ts.d
|
|---|
| 349 | do i=1,num_rep
|
|---|
| 350 | j=intem(i)
|
|---|
| 351 | write(14,*) iold, i, j, pbe(i),
|
|---|
| 352 | & yol(j), eyslr(j), eysmip(j),
|
|---|
| 353 | & rgyrp(j), nhelp(j), nbetp(j), imhbp(j)
|
|---|
| 354 |
|
|---|
| 355 | end do
|
|---|
| 356 | c Write the current parallel tempering information into par_R.in
|
|---|
| 357 | if ((mod(iold, nsave).eq.0))
|
|---|
| 358 | & then
|
|---|
| 359 | open(13,file='par_R.in', status='unknown')
|
|---|
| 360 | write(13,*) iold
|
|---|
| 361 | do i=1,num_rep
|
|---|
| 362 | write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),
|
|---|
| 363 | & h_maxp(i)
|
|---|
| 364 | end do
|
|---|
| 365 | C -------------------------- Various statistics of current run
|
|---|
| 366 | c swp=nswp-nequi
|
|---|
| 367 | swp=nsw
|
|---|
| 368 | write(13,*) 'Acceptance rate for change of chains:'
|
|---|
| 369 | do k1=1,num_rep
|
|---|
| 370 | temp=1.0d0/pbe(k1)/0.00198773
|
|---|
| 371 | write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
|
|---|
| 372 | c Above: it's the acceptance rate of exchange of replicas. Since a
|
|---|
| 373 | c replica exchange is attempted only once every nmes sweeps, the
|
|---|
| 374 | c rate should be normalized with (nmes/swp).
|
|---|
| 375 | end do
|
|---|
| 376 | write(13,*)
|
|---|
| 377 | do k1=1,num_rep
|
|---|
| 378 | k = intem(k1)
|
|---|
| 379 | temp=1.0d0/pbe(k1)/0.00198773
|
|---|
| 380 | beta = pbe(k1)
|
|---|
| 381 | geavm(k1) = nmes*eavm(k1)/swp
|
|---|
| 382 | gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2)
|
|---|
| 383 | # *beta*beta/nresi
|
|---|
| 384 | write(13,'(a,2f9.2,i4,f12.3)')
|
|---|
| 385 | & 'Temperature, Node,local acceptance rate:',
|
|---|
| 386 | & beta,temp,k,acy(k1)/dble(nsw*nvr)
|
|---|
| 387 | c Above: Changed (nswp-nequi) in the denominator of acceptance as
|
|---|
| 388 | c acceptance values are initialized to 0 after equilibration cycles are
|
|---|
| 389 | c finished. Note also that since this is being written in the middle of
|
|---|
| 390 | c the simulation, it is normalized to nsw instead of nswp.
|
|---|
| 391 | write(13,'(a,3f12.2)')
|
|---|
| 392 | & 'Last Energy, Average Energy, Spec. Heat:',
|
|---|
| 393 | & yol(k),geavm(k1),gsph(k1)
|
|---|
| 394 | write(13,*)
|
|---|
| 395 | end do
|
|---|
| 396 | close(13)
|
|---|
| 397 | ! Finally, flush the time series file to ensure that we can do a proper
|
|---|
| 398 | ! restart.
|
|---|
| 399 | call flush(14)
|
|---|
| 400 | call flush(18)
|
|---|
| 401 | end if
|
|---|
| 402 |
|
|---|
| 403 | C--------------------Parallel Tempering update
|
|---|
| 404 | c Swap with right neighbor (odd, even)
|
|---|
| 405 | if(odd.eq.1) then
|
|---|
| 406 | nu=1
|
|---|
| 407 | no1 = num_rep-1
|
|---|
| 408 | c Swap with left neighbor (even, odd)
|
|---|
| 409 | else
|
|---|
| 410 | nu = 2
|
|---|
| 411 | no1 = num_rep
|
|---|
| 412 | end if
|
|---|
| 413 | do i=nu,no1,2
|
|---|
| 414 | j=i+1
|
|---|
| 415 | c Periodic bc for swaps
|
|---|
| 416 | if(i.eq.num_rep) j=1
|
|---|
| 417 | iidx=intem(i)
|
|---|
| 418 | jn=intem(j)
|
|---|
| 419 | wij=exp(-pbe(i)*yol(jn)-pbe(j)*yol(iidx)
|
|---|
| 420 | & +pbe(i)*yol(iidx)+pbe(j)*yol(jn))
|
|---|
| 421 | rd=grnd()
|
|---|
| 422 | randomCount = randomCount + 1
|
|---|
| 423 | if(wij.ge.rd) then
|
|---|
| 424 | acx1(i) = acx1(i)+1
|
|---|
| 425 | intem(i) = jn
|
|---|
| 426 | intem(j) = iidx
|
|---|
| 427 | inode(in)= j
|
|---|
| 428 | inode(jn)= i
|
|---|
| 429 | end if
|
|---|
| 430 | end do
|
|---|
| 431 | c ---------------- End Loop over nodes which creates a new temperature
|
|---|
| 432 | c map for all nodes, at the node with rank 0.
|
|---|
| 433 | c
|
|---|
| 434 | odd = 1 - odd
|
|---|
| 435 | end if
|
|---|
| 436 | c End of "if (myrank.eq.0) ...". The block above includes PT update and
|
|---|
| 437 | c writing of observables into the time series file etc.
|
|---|
| 438 |
|
|---|
| 439 | c Below: Communicate new temperature-node map to all nodes
|
|---|
| 440 | CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
|
|---|
| 441 | & IERR)
|
|---|
| 442 | CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
|
|---|
| 443 | & IERR)
|
|---|
| 444 | CALL MPI_BCAST(E_MINP,num_rep,MPI_DOUBLE_PRECISION,0,
|
|---|
| 445 | & MPI_COMM_WORLD,IERR)
|
|---|
| 446 | CALL MPI_BCAST(H_MAXP,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
|
|---|
| 447 | & IERR)
|
|---|
| 448 | c Synchronize random number generators for replica 0
|
|---|
| 449 | if (rep_id.eq.0) then
|
|---|
| 450 | CALL MPI_BCAST(randomCount,1,MPI_INTEGER,0,my_mpi_comm,
|
|---|
| 451 | & IERR)
|
|---|
| 452 | if (myrank.ne.0) then
|
|---|
| 453 | ! write (*,*) '[', myrank,'] Missed', randomCount,
|
|---|
| 454 | ! & 'random numbers.'
|
|---|
| 455 | do i = 1, randomCount
|
|---|
| 456 | rd = grnd()
|
|---|
| 457 | ! write (*,*) '[', myrank,'] rd=', rd
|
|---|
| 458 | enddo
|
|---|
| 459 | endif
|
|---|
| 460 | endif
|
|---|
| 461 |
|
|---|
| 462 | BETA=PBE(INODE(rep_id+1))
|
|---|
| 463 | e_min = e_minp(inode(rep_id+1))
|
|---|
| 464 | h_max = h_maxp(inode(rep_id+1))
|
|---|
| 465 | if (INODE(rep_id + 1).eq.1) dir = 1
|
|---|
| 466 | if (INODE(rep_id + 1).eq.num_rep) dir = -1
|
|---|
| 467 |
|
|---|
| 468 | endif
|
|---|
| 469 | c End of "if (mod(iold,nmes).eq.0) ..."
|
|---|
| 470 | end do
|
|---|
| 471 | c-----------End Loop over sweeps
|
|---|
| 472 | c
|
|---|
| 473 | C OUTPUT:
|
|---|
| 474 | C--------------------For Re-starts:
|
|---|
| 475 | nu = rep_id + 1
|
|---|
| 476 | filebase = "conf_0000.var"
|
|---|
| 477 | call outvar(0, fileNameMP(filebase, 6, 9, nu))
|
|---|
| 478 | e_final=energy()
|
|---|
| 479 | if (partem_comm.ne.MPI_COMM_NULL) then
|
|---|
| 480 | write (*,*) rep_id, ' E_final', e_final
|
|---|
| 481 | endif
|
|---|
| 482 | eol0 = eol
|
|---|
| 483 | acz0 = acz
|
|---|
| 484 | if (partem_comm.ne.MPI_COMM_NULL) then
|
|---|
| 485 | CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,
|
|---|
| 486 | # MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 487 | CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1,
|
|---|
| 488 | # MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
|
|---|
| 489 | endif
|
|---|
| 490 |
|
|---|
| 491 | if(rep_id.eq.0.and.myrank.eq.0) then
|
|---|
| 492 | close(14)
|
|---|
| 493 | open(13,file='par_R.in', status='unknown')
|
|---|
| 494 | write(13,*) iold
|
|---|
| 495 | do i=1,num_rep
|
|---|
| 496 | write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
|
|---|
| 497 | end do
|
|---|
| 498 | C -------------------------- Various statistics of current run
|
|---|
| 499 | swp=nswp
|
|---|
| 500 | write(13,*) 'Acceptance rate for change of chains:'
|
|---|
| 501 | do k1=1,num_rep
|
|---|
| 502 | temp=1.0d0/pbe(k1)/0.00198773
|
|---|
| 503 | write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
|
|---|
| 504 | end do
|
|---|
| 505 | write(13,*)
|
|---|
| 506 | do k1=1,num_rep
|
|---|
| 507 | k = intem(k1)
|
|---|
| 508 | temp=1.0d0/pbe(k1)/0.00198773
|
|---|
| 509 | beta = pbe(k1)
|
|---|
| 510 | geavm(k1) = nmes*eavm(k1)/swp
|
|---|
| 511 | gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2)*beta*beta/nresi
|
|---|
| 512 | write(13,'(a,2f9.2,i4,f12.3)')
|
|---|
| 513 | & 'Temperature, Node,local acceptance rate:',
|
|---|
| 514 | & beta,temp,k,acy(k1)/dble((nswp)*nvr)
|
|---|
| 515 | write(13,'(a,3f12.2)')
|
|---|
| 516 | & 'Last Energy, Average Energy, Spec. Heat:',
|
|---|
| 517 | & yol(k),geavm(k1),gsph(k1)
|
|---|
| 518 | write(13,*)
|
|---|
| 519 | end do
|
|---|
| 520 | close(13)
|
|---|
| 521 | c close(16)
|
|---|
| 522 | end if
|
|---|
| 523 | close(18)
|
|---|
| 524 |
|
|---|
| 525 | c =====================
|
|---|
| 526 | CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
|
|---|
| 527 |
|
|---|
| 528 | return
|
|---|
| 529 |
|
|---|
| 530 | end
|
|---|
