Changeset bd2278d for EXAMPLES/partem_p.f
- Timestamp:
- 09/05/08 11:49:42 (16 years ago)
- Branches:
- master
- Children:
- fafe4d6
- Parents:
- 2ebb8b6
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
EXAMPLES/partem_p.f
r2ebb8b6 rbd2278d 1 c**************************************************************2 c3 cThis file contains the subroutines: partem_p4 CCompared to the version in the main distribution, this5 Croutine doesn't write the rmsd nor native contacts to the time6 Cseries.7 c8 cCopyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,9 cShura Hayryan, Chin-Ku Hu10 cCopyright 2007 Frank Eisenmenger, U.H.E. Hansmann,11 cJan H. Meinke, Sandipan Mohanty12 c13 c**************************************************************1 !************************************************************** 2 ! 3 ! This file contains the subroutines: partem_p 4 ! Compared to the version in the main distribution, this 5 ! routine doesn't write the rmsd nor native contacts to the time 6 ! series. 7 ! 8 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 9 ! Shura Hayryan, Chin-Ku Hu 10 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 11 ! Jan H. Meinke, Sandipan Mohanty 12 ! 13 ! ************************************************************** 14 14 15 15 subroutine partem_p(num_rep, nequi, nswp, nmes, nsave, newsta, 16 16 & switch, rep_id, partem_comm) 17 C18 CPURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM19 CON PARALLEL COMPUTERS USING MPI20 C21 Cswitch: Choses the starting configuration:22 C-1 - stretched configuration23 C0 - don't change anything24 C1 - random start configuration25 C26 cCALLS: addang,contacts,energy,hbond,helix,iendst,metropolis,27 coutvar,(rand),rgyr28 C17 ! 18 ! PURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM 19 ! ON PARALLEL COMPUTERS USING MPI 20 ! 21 ! switch: Choses the starting configuration: 22 ! -1 - stretched configuration 23 ! 0 - don't change anything 24 ! 1 - random start configuration 25 ! 26 ! CALLS: addang,contacts,energy,hbond,helix,iendst,metropolis, 27 ! outvar,(rand),rgyr 28 ! 29 29 include '../INCL.H' 30 30 include '../INCP.H' … … 37 37 external can_weight 38 38 39 Cnequi: number of Monte Carlo sweeps for thermalization40 Cnswp: number of Monte Carlo sweeps41 Cnmes: number of Monte Carlo sweeps between measurments42 Cnewsta: .true. for new simulations, .false. for re-start39 ! nequi: number of Monte Carlo sweeps for thermalization 40 ! nswp: number of Monte Carlo sweeps 41 ! nmes: number of Monte Carlo sweeps between measurments 42 ! newsta: .true. for new simulations, .false. for re-start 43 43 44 44 dimension eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC), … … 50 50 double precision e_min, e_minp(MAX_PROC), e_minpt(MAX_PROC) 51 51 integer h_max, h_maxp(MAX_PROC) 52 cOrder of replica exchange52 ! Order of replica exchange 53 53 integer odd 54 54 ! Counter to keep random number generators in sync 55 55 integer randomCount 56 56 57 cCollect partial energies. Only the root writes to disk. We have to58 ccollect the information from the different replicas and provide59 carrays to store them.60 ceyslr storage array for solvent energy61 ceyelp - " - coulomb energy62 ceyvwp - " - van-der-Waals energy63 ceyhbp - " - hydrogen bonding energy64 ceysmi - " - intermolecular interaction energy57 ! Collect partial energies. Only the root writes to disk. We have to 58 ! collect the information from the different replicas and provide 59 ! arrays to store them. 60 ! eyslr storage array for solvent energy 61 ! eyelp - " - coulomb energy 62 ! eyvwp - " - van-der-Waals energy 63 ! eyhbp - " - hydrogen bonding energy 64 ! eysmi - " - intermolecular interaction energy 65 65 double precision eyslr(MAX_PROC) 66 66 double precision eyelp(MAX_PROC),eyvwp(MAX_PROC),eyhbp(MAX_PROC), 67 67 & eyvrp(MAX_PROC),eysmip(MAX_PROC) 68 cCollect information about accessible surface and van-der-Waals volume69 casap storage array for solvent accessible surface70 cvdvolp storage array for van-der-Waals volume68 ! Collect information about accessible surface and van-der-Waals volume 69 ! asap storage array for solvent accessible surface 70 ! vdvolp storage array for van-der-Waals volume 71 71 double precision asap(MAX_PROC), vdvolp(MAX_PROC) 72 72 … … 75 75 integer imhbp(MAX_PROC) 76 76 character*80 filebase, fileNameMP, tbase0,tbase1 77 cframe frame number for writing configurations78 ctrackID configuration that should be tracked and written out79 cdir direction in random walk80 c-1 - visited highest temperature last81 c1 - visited lowest temperature last82 c0 - haven't visited the boundaries yet.83 cdirp storage array for directions.77 ! frame frame number for writing configurations 78 ! trackID configuration that should be tracked and written out 79 ! dir direction in random walk 80 ! -1 - visited highest temperature last 81 ! 1 - visited lowest temperature last 82 ! 0 - haven't visited the boundaries yet. 83 ! dirp storage array for directions. 84 84 integer frame, trackID, dir 85 85 integer dirp(MAX_PROC) … … 92 92 & rep_id, num_rep, partem_comm, myrank 93 93 call flush(6) 94 C95 c96 CFile with temperatures94 ! 95 ! 96 ! File with temperatures 97 97 open(11,file='temperatures_abeta',status='old') 98 98 … … 100 100 open(18,file=fileNameMP(tbase0,5,9,rep_id),status='unknown') 101 101 if (rep_id.eq.0.and.myrank.eq.0) then 102 cFile with time series of simulation102 ! File with time series of simulation 103 103 open(14,file='ts.d',status='unknown') 104 104 endif 105 105 106 CREAD IN TEMPERATURES106 ! READ IN TEMPERATURES 107 107 do i=1,num_rep 108 108 read(11,*) j,temp … … 111 111 close(11) 112 112 113 cnresi: number of residues113 ! nresi: number of residues 114 114 nresi=irsml2(1)-irsml1(1)+1 115 C116 CInitialize variables115 ! 116 ! Initialize variables 117 117 do i=1,num_rep 118 118 acx1(i) = 0.0d0 … … 132 132 dir = dirp(rep_id + 1) 133 133 134 c_________________________________ Initialize Variables134 ! _________________________________ Initialize Variables 135 135 if(newsta) then 136 136 iold=0 … … 139 139 intem(i) = i 140 140 end do 141 c_________________________________ initialize starting configuration141 ! _________________________________ initialize starting configuration 142 142 if (switch.ne.0) then 143 143 do i=1,nvr … … 173 173 CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,IERR) 174 174 CALL MPI_BCAST(YOL,num_rep,MPI_DOUBLE_PRECISION,0, 175 #MPI_COMM_WORLD,IERR)175 & MPI_COMM_WORLD,IERR) 176 176 CALL MPI_BCAST(E_MINP, num_rep, MPI_DOUBLE_PRECISION, 0, 177 #MPI_COMM_WORLD, IERR)177 & MPI_COMM_WORLD, IERR) 178 178 CALL MPI_BCAST(h_maxp,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD, 179 $IERR)179 & IERR) 180 180 end if 181 181 … … 189 189 write(*,*) rep_id, yol(rep_id + 1), eol 190 190 endif 191 CStart of simulation191 ! Start of simulation 192 192 write (*,*) '[',rep_id, myrank, beta, partem_comm, 193 193 & '] Energy before equilibration:', eol 194 c=====================Equilibration by canonical Metropolis194 ! =====================Equilibration by canonical Metropolis 195 195 do nsw=1,nequi 196 196 call metropolis(eol,acz,can_weight) … … 199 199 write (*,*) '[',rep_id,'] Energy after equilibration:', eol 200 200 call flush(6) 201 C202 C======================Multiple Markov Chains201 ! 202 !======================Multiple Markov Chains 203 203 acz = 0 204 204 do nsw=1,nswp 205 c------------First ordinary Metropolis205 !------------First ordinary Metropolis 206 206 call metropolis(eol,acz,can_weight) 207 207 iold = iold + 1 … … 214 214 endif 215 215 acz0 = acz 216 CMeasure global radius of gyration216 ! Measure global radius of gyration 217 217 call rgyr(0,rgy,ee) 218 218 rgyp = rgy 219 CMeasure Helicity and Sheetness219 ! Measure Helicity and Sheetness 220 220 call helix(nhel,mhel,nbet,mbet) 221 CMeasure Number of hydrogen bonds221 ! Measure Number of hydrogen bonds 222 222 mhb = 0 223 223 do i = 1, ntlml … … 226 226 enddo 227 227 call interhbond(imhb) 228 CMeasure total number of contacts (NCTOT) and number of229 Cnative contacts (NCNAT)228 ! Measure total number of contacts (NCTOT) and number of 229 ! native contacts (NCNAT) 230 230 call contacts(nctot,ncnat,dham) 231 cAdd tracking of lowest energy configuration231 ! Add tracking of lowest energy configuration 232 232 if (eol.lt.e_min) then 233 cWrite out configuration233 ! Write out configuration 234 234 i=rep_id+1 235 235 j=inode(i) … … 248 248 close(15) 249 249 endif 250 cAdd tracking of configuration with larges hydrogen contents.250 ! Add tracking of configuration with larges hydrogen contents. 251 251 if ((mhb + imhb).gt.h_max) then 252 cWrite out configuration252 ! Write out configuration 253 253 i = rep_id + 1 254 254 j = inode(i) … … 268 268 endif 269 269 270 C271 C--------------------Gather measurement data270 ! 271 !--------------------Gather measurement data 272 272 ! I only use the master node of each replica for data collection. The 273 273 ! variable partem_comm provides the appropriate communicator. … … 310 310 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR) 311 311 312 cWrite trajectory312 ! Write trajectory 313 313 write (18,*) '@@@',iold,inode(rep_id+1) 314 314 call outvbs(0,18) 315 315 write (18,*) '###' 316 316 ! call flush(18) 317 cWrite current configuration317 ! Write current configuration 318 318 if ((mod(iold, nsave).eq.0)) then 319 319 filebase = "conf_0000.var" … … 324 324 if(rep_id.eq.0.and.myrank.eq.0) then 325 325 randomCount = 0 326 cUpdate acceptance, temperature wise average of E and E^2 used to calculate327 cspecific heat.326 ! Update acceptance, temperature wise average of E and E^2 used to calculate 327 ! specific heat. 328 328 do i=1,num_rep 329 329 j=intem(i) 330 330 acy(i)=0.0 331 cAbove: contents of acy1 are added to acy(i) a few lines down.332 cacy1(intem(i)) contains information received from the node at temperature333 ci, on how many updates have been accepted in node intem(i). Since acz334 cis not reset to 0 every cycle, acy(i) must be set to 0 here. Else, there335 cwill be serious double counting and the values of acceptance printed336 cwill be simply wrong.331 ! Above: contents of acy1 are added to acy(i) a few lines down. 332 ! acy1(intem(i)) contains information received from the node at temperature 333 ! i, on how many updates have been accepted in node intem(i). Since acz 334 ! is not reset to 0 every cycle, acy(i) must be set to 0 here. Else, there 335 ! will be serious double counting and the values of acceptance printed 336 ! will be simply wrong. 337 337 e_minpt(i)=e_minp(intem(i)) 338 338 end do … … 346 346 347 347 348 CWrite measurements to the time series file ts.d348 ! Write measurements to the time series file ts.d 349 349 do i=1,num_rep 350 350 j=intem(i) … … 354 354 355 355 end do 356 cWrite the current parallel tempering information into par_R.in356 ! Write the current parallel tempering information into par_R.in 357 357 if ((mod(iold, nsave).eq.0)) 358 358 & then … … 363 363 & h_maxp(i) 364 364 end do 365 C-------------------------- Various statistics of current run366 cswp=nswp-nequi365 ! -------------------------- Various statistics of current run 366 ! swp=nswp-nequi 367 367 swp=nsw 368 368 write(13,*) 'Acceptance rate for change of chains:' … … 370 370 temp=1.0d0/pbe(k1)/0.00198773 371 371 write(13,*) temp, acx1(k1)*2.0d0*nmes/swp 372 cAbove: it's the acceptance rate of exchange of replicas. Since a373 creplica exchange is attempted only once every nmes sweeps, the374 crate should be normalized with (nmes/swp).372 ! Above: it's the acceptance rate of exchange of replicas. Since a 373 ! replica exchange is attempted only once every nmes sweeps, the 374 ! rate should be normalized with (nmes/swp). 375 375 end do 376 376 write(13,*) … … 381 381 geavm(k1) = nmes*eavm(k1)/swp 382 382 gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2) 383 #*beta*beta/nresi383 & *beta*beta/nresi 384 384 write(13,'(a,2f9.2,i4,f12.3)') 385 385 & 'Temperature, Node,local acceptance rate:', 386 386 & beta,temp,k,acy(k1)/dble(nsw*nvr) 387 cAbove: Changed (nswp-nequi) in the denominator of acceptance as388 cacceptance values are initialized to 0 after equilibration cycles are389 cfinished. Note also that since this is being written in the middle of390 cthe simulation, it is normalized to nsw instead of nswp.387 ! Above: Changed (nswp-nequi) in the denominator of acceptance as 388 ! acceptance values are initialized to 0 after equilibration cycles are 389 ! finished. Note also that since this is being written in the middle of 390 ! the simulation, it is normalized to nsw instead of nswp. 391 391 write(13,'(a,3f12.2)') 392 392 & 'Last Energy, Average Energy, Spec. Heat:', … … 401 401 end if 402 402 403 C--------------------Parallel Tempering update404 cSwap with right neighbor (odd, even)403 !--------------------Parallel Tempering update 404 ! Swap with right neighbor (odd, even) 405 405 if(odd.eq.1) then 406 406 nu=1 407 407 no1 = num_rep-1 408 cSwap with left neighbor (even, odd)408 ! Swap with left neighbor (even, odd) 409 409 else 410 410 nu = 2 … … 413 413 do i=nu,no1,2 414 414 j=i+1 415 cPeriodic bc for swaps415 ! Periodic bc for swaps 416 416 if(i.eq.num_rep) j=1 417 417 iidx=intem(i) … … 429 429 end if 430 430 end do 431 c---------------- End Loop over nodes which creates a new temperature432 cmap for all nodes, at the node with rank 0.433 c431 ! ---------------- End Loop over nodes which creates a new temperature 432 ! map for all nodes, at the node with rank 0. 433 ! 434 434 odd = 1 - odd 435 435 end if 436 cEnd of "if (myrank.eq.0) ...". The block above includes PT update and437 cwriting of observables into the time series file etc.436 ! End of "if (myrank.eq.0) ...". The block above includes PT update and 437 ! writing of observables into the time series file etc. 438 438 439 cBelow: Communicate new temperature-node map to all nodes439 ! Below: Communicate new temperature-node map to all nodes 440 440 CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD, 441 441 & IERR) … … 446 446 CALL MPI_BCAST(H_MAXP,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD, 447 447 & IERR) 448 cSynchronize random number generators for replica 0448 ! Synchronize random number generators for replica 0 449 449 if (rep_id.eq.0) then 450 450 CALL MPI_BCAST(randomCount,1,MPI_INTEGER,0,my_mpi_comm, … … 467 467 468 468 endif 469 cEnd of "if (mod(iold,nmes).eq.0) ..."469 ! End of "if (mod(iold,nmes).eq.0) ..." 470 470 end do 471 c-----------End Loop over sweeps472 c473 COUTPUT:474 C--------------------For Re-starts:471 !-----------End Loop over sweeps 472 ! 473 ! OUTPUT: 474 !--------------------For Re-starts: 475 475 nu = rep_id + 1 476 476 filebase = "conf_0000.var" … … 484 484 if (partem_comm.ne.MPI_COMM_NULL) then 485 485 CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1, 486 #MPI_DOUBLE_PRECISION,0,partem_comm,IERR)486 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR) 487 487 CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1, 488 #MPI_DOUBLE_PRECISION,0,partem_comm,IERR)488 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR) 489 489 endif 490 490 … … 496 496 write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i) 497 497 end do 498 C-------------------------- Various statistics of current run498 ! -------------------------- Various statistics of current run 499 499 swp=nswp 500 500 write(13,*) 'Acceptance rate for change of chains:' … … 519 519 end do 520 520 close(13) 521 cclose(16)521 ! close(16) 522 522 end if 523 523 close(18) 524 524 525 c=====================525 ! ===================== 526 526 CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) 527 527
Note:
See TracChangeset
for help on using the changeset viewer.