1* MDynaMix v. > 5.1 2* 3* File pimd.f 4* ------------ 5* 6C This file contains subroutines relevant to Path Integral Molecular Dynamics 7* 8*==================================================================== 9* 10* 1. PIMD - Integration of motion in Path Integral Molecular Dynamics 11* 12* Centroid MD in the adiabatic approximation and 13* use of double time step algorithm 14* (will be described in more details in some future) 15* 16* This subroutine repeats in many features the DOUBLE subroutine in mdstep.f file 17C 18 SUBROUTINE PIMD 19 include "prcm.h" 20* 21* 1.1 Local definitions 22* 23C LRDF defines wheter to calculate RDFs 24 LOGICAL LRDF,DONE 25C this is temporary array to keep long-range contributions to virial 26 real*8 TRLL(3) 27C RDF should not be calculated during initialisation stage (2.2) 28* write(*,*)' node ',JBEAD,' enter PIMD' 29 LRDF=.false. 30 NFREQ2=NFREQ/2+1 31 VIRA(1)=0. 32 VIRA(2)=0. 33 VIRA(3)=0. 34 NEKS=0 35 QEKS=0. 36 TEMPNOD=0. 37 IPIMA=0 38 do I=1,NSITES 39 PIRS(I)=0. 40 end do 41 do I=1,NSTOT 42 do J=1,JNH 43 THERMS(I,J)=0. 44 end do 45 end do 46* 47* 1.2 Preparation 48* --------------- 49* 50* 1.2.1 Compose list of neigbours 51 call CHCNB(LRDF) ! mdstep.f 52* 1.2.2 Calculate centre-of-masses 53 CALL GETCOM 54* 1.2.3 Calculate all forces (should be done before entering the cycle) 55 LRDF=.true. 56C slow forces (important: slow forces calculated before fast ones!) 57 call SLFORCES 58C Fast forces 59 call FFORCES(.false.) 60C Scale down the forces 61 do I=1,NSTOT 62 GX(I)=GX(I)/NBEADS 63 GY(I)=GY(I)/NBEADS 64 GZ(I)=GZ(I)/NBEADS 65 HX(I)=HX(I)/NBEADS 66 HY(I)=HY(I)/NBEADS 67 HZ(I)=HZ(I)/NBEADS 68 TEMBS(I)=0. 69 end do 70C Forces between PI beads 71 call PIFORCES 72* write(*,'(i3,a8,7e12.5)')JBEAD,' PIRR = ',(PIRR(IS),IS=1,NSITES) 73* 1.2.4 Calculate different contributions to pressure and energy 74C (this is done mostly for control purposes) 75 timeg0 = cputime(0.d0) 76 if(IPRINT.ge.7)then 77 PELL=PELS1+PELS2+SPE 78 PENL=(PE1+PE2+PELL)*PERMOL 79 PEL=PELL*PERMOL 80 PENS=PINT*PERMOL 81 PEKI=TKE*PERMOL 82 ETOT=PENL+PENS+PEKI 83 WRITE(6,'(4(a6,f12.5))') 84 +' Eex= ',PENL,' Ein= ',PENS,' Eel= ',PEL,' Etot=',ETOT 85 end if 86* 87* 1.3 Double time step algorithm 88* ------------------------------- 89* 1.3.1 Organize the cycle (till 2.4.9) 90 NSTEP=0 91 1 NSTEP=NSTEP+1 92 MSTEP=NSTTOT+NSTEP 93* 1.3.2 First thermostat integration: to half of the long time step 94 call PI_SCALING(1) 95* 1.3.3 First integration of slow forces 96C Velocities are corrected as due to slow forces acting during 97C half of the long time step 98 DO I = NAB(TASKID),NAE(TASKID) 99 ITYP = ITYPE(I) 100C LMOVE define which molecules can move (normally, LMOVE=.true.) 101 if(LMOVE(ITYP))then 102 VX(I) = VX(I)+HTSTEP*GX(I) 103 VY(I) = VY(I)+HTSTEP*GY(I) 104 VZ(I) = VZ(I)+HTSTEP*GZ(I) 105 end if 106 END DO! OF I 107* 108* 1.3.4 Integration of the fast forces 109* ------------------------------------- 110* 1.3.4.1 Cycle over short time steps 111C This set zeros for current averages bond/angles and their energies 112C These averages are taken over the molecules and over short time 113C steps within this large step 114 call ZEROAV 115 QEKS=0 116 DO NSTEB = 1,NFREQ 117 DO I = NAB(TASKID),NAE(TASKID) 118 ITYP = ITYPE(I) 119 if(LMOVE(ITYP))then 120* 1.3.4.2 Velocities at t(short)+1/2 121 VX(I) = VX(I)+HTSTEB*HX(I) ! V(t+1/2) 122 VY(I) = VY(I)+HTSTEB*HY(I) 123 VZ(I) = VZ(I)+HTSTEB*HZ(I) 124 end if 125 end do 126* 1.3.4.3. Computation of centroid momentum FX,FY,FZ 127* call CMD_MOM ! mpi communication 128* 1.3.4.4 Nose themostat for local motions of beads 129 call BEADNOSE 130* 1.3.4.5 Coordinates at t(short)+1 131 DO I = NAB(TASKID),NAE(TASKID) 132 ITYP = ITYPE(I) 133 if(LMOVE(ITYP))then 134 SX(I) = VX(I)*TSTEB*MASSDI(I)+SX(I) ! X(t+1) 135 SY(I) = VY(I)*TSTEB*MASSDI(I)+SY(I) 136 SZ(I) = VZ(I)*TSTEB*MASSDI(I)+SZ(I) 137 end if 138 END DO! OF I 139 timeg = timeg+cputime(timeg0) 140* 1.3.4.6 Recalculate centre-of-mass coordinates 141 CALL GETCOM ! services.f 142* 143* 1.3.5 Recalculate forces 144* ------------------------ 145* 1.3.5.1 After "ICHNB" long time steps - recalculate list of neighbours 146C This is done in the middle of long time step 147 if(NSTEB.eq.NFREQ2.and.mod(MSTEP,ICHNB).eq.0)call CHCNB(LRDF) 148 if(NSTEB.eq.NFREQ)then 149* 1.3.5.2 Each long time step - calculate slow forces 150 CALL SLFORCES 151* 1.3.5.3 Claculate fast forces 152C The logical parameter signals whether to calculate averages 153 call FFORCES(.true.) 154 do I=1,NSTOT 155 GX(I)=GX(I)/NBEADS 156 GY(I)=GY(I)/NBEADS 157 GZ(I)=GZ(I)/NBEADS 158 HX(I)=HX(I)/NBEADS 159 HY(I)=HY(I)/NBEADS 160 HZ(I)=HZ(I)/NBEADS 161 end do 162C Forces between PI beads 163 call PIFORCES 164 else 165 call FFORCES(.false.) 166 do I=1,NSTOT 167 HX(I)=HX(I)/NBEADS 168 HY(I)=HY(I)/NBEADS 169 HZ(I)=HZ(I)/NBEADS 170 end do 171 call PIFORCES 172 end if 173 QEKS = QEKS + QEK 174* 175* 1.3.6 Conclude integration of fast forces 176* ----------------------------------------- 177C (continued in the beginning of the cycle, see p.2.3.2 178 timeg0 = cputime(0.d0) 179 if(IPRINT.ge.8)write(*,'(a,I3,f12.3)')' 1.st.',NSTEB,TEMP 180 DO I = NAB(TASKID),NAE(TASKID) 181 ITYP = ITYPE(I) 182 if(LMOVE(ITYP))then 183 VX(I) = VX(I)+HTSTEB*HX(I) ! V(t+1) 184 VY(I) = VY(I)+HTSTEB*HY(I) 185 VZ(I) = VZ(I)+HTSTEB*HZ(I) 186 end if 187 END DO! OF I 188 call GETEMP 189 if(IPRINT.ge.8)write(*,'(a,2I3,f12.3)')' sh.st.',JBEAD,NSTEB,TEMP 190* 191 END DO! OF NSTEB 192* 193* 1.3.7 Conclude integration of slow forces 194* ----------------------------------------- 195* 1.3.7.1 Velocities at t(long)+1 196 DO I = NAB(TASKID),NAE(TASKID) 197 ITYP = ITYPE(I) 198 if(LMOVE(ITYP))then 199 VX(I) = VX(I)+HTSTEP*GX(I) 200 VY(I) = VY(I)+HTSTEP*GY(I) 201 VZ(I) = VZ(I)+HTSTEP*GZ(I) 202 end if 203 END DO! OF I 204 timeg = timeg+cputime(timeg0) 205* 1.3.7.2 Conclude integration of NHC 206 call PI_SCALING(2) 207* 1.3.6.3 Report velocities to the master node 208 timeg0 = cputime(0.d0) 209* 210* 1.4 calculate averages for intermediate output 211* ---------------------------------------------- 212* remove excess COM momenta if specified 213 if(ICMOM.gt.0.and.mod(MSTEP,ICMOM).eq.0)call CHKMOM 214* 1.4.1 Different contribution to the kinetic energy 215 CALL GETKIN ! services.f 216* Accumullation of data from the nodes 217* 1.4.4 Acccumulate averages 218 TEMPNOD=TEMPNOD+TEMP 219 call PI_AVER 220 PEST = PELS1+PELS2+SPE 221 PE = PE1+PE2+PEST 222 if(mod(MSTEP,IAVER).eq.0)CALL GETAVR(1) ! aver.f 223* 1.4.5 Dump trajectory 224 if(mod(MSTEP,NTREK).eq.0)CALL TRACE ! restart.f 225* 1.4.6 Calculate different contributions to pressure and energy 226 PENL = PE*PERMOL ! intermolecular energy 227 PEL= PEST*PERMOL ! electrostatic energy 228 PENS=PINT*PERMOL ! intramolecular energy 229 PEKI=TKE*PERMOL ! classical kinetic energy 230 POTEN=PENL+PENS ! total potential energy 231 ETOT =POTEN+PEKI ! total energy 232 WIRSUM=WIRS+WIRSS 233 TRYCKE = PEL*UNITP/(3.*VOL*PERMOL) ! electrostatic 234 TRYCKL = (VIR1+VIR2)*UNITP/(3.*VOL) ! LJ 235 TRYCKB = VIRB*UNITP/(3.*VOL) ! bonds 236 TRYCKS = WIRSUM*UNITP/(3.*VOL) ! molecular 237 TRYCKP = VIRPI*UNITP/(3.*VOL) ! PI springs 238* TRYCKA = (VIRA(1)+VIRA(2)+VIRA(3))*UNITP/(3.*VOL) ! molecular 239 TRYCKD = VIRD*UNITP/(3.*VOL) ! long-range corr. 240 DENS = TOTMAS*1.d-3/(VOL*UNITL**3) 241C Quantum kinetic energy 242 QEK = QEKS / NFREQ 243 TKEQ = 1.5*NSTOT*NBEADS*BOLTZ*TRTEMP/UNITE-QEK 244* 1.4.7 Current output 245C OUTPUT for visualisation 246 if(LVISUAL)then 247 if(mod(MSTEP,IOUTPUT).eq.0)then 248 write(*,'(a)')'@mm coord.start' 249 do I=1,NSTOT 250 write(*,'(a,3f12.4)')'@mm ',SX(I),SY(I),SZ(I) 251 end do 252 write(*,'(a)')'@mm coord.end' 253 end if 254 WRITE(*,'(a,I7,4(1X,F9.3),2(1X,f9.2),f9.4,3f10.4)') 255 + '@mm E ',MSTEP,POTEN,PENS,PEL,ETOT,TEMP,TRYCKM,DENS, 256 + BOXL,BOYL,BOZL 257 end if ! if(LVISUAL 258 if(mod(MSTEP,IOUTPUT).eq.0)then 259 if(IPRINT.ge.5)then 260 WRITE(6,'(I8,4(1X,F9.3),2(1X,f9.2),f9.4)') 261 + MSTEP,POTEN,PENS,PEL,ETOT,TEMP,TRYCK,DENS 262 if(JBEAD.eq.0)write(*,'(a,f10.3,a,f10.3)') 263 + 'Eharm: ',QEK*PERMOL,' Ekin_Q:',TKEQ*PERMOL 264 end if 265 if(IPRINT.ge.6)then 266 write(*,'(2f7.4,f9.2,4f7.0,I7,I9)') 267 + SC,SCL,TRYCKM,TTR,TROT,TINT,TEMPV,MBSH,MBLN 268 if(JBEAD.eq.0)write(*,'(a,80f8.4)') 269 + 'RRav: ',(PIRR(I),I=1,NSITES) 270 if(LSEP)write(*,'(4f11.1,3f11.3)') 271 + TRYCKX,TRYCKY,TRYCKZ,(TRYCKX+TRYCKY+TRYCKZ)/3.,BOXL,BOYL,BOZL 272 if(IEXT.eq.1)write(*,'(2(a,e14.6))') 273 + ' E abs: ',(EABS+EABSS)*PERMOL,' Ext.f. ',EFIELD 274 end if 275 if(IPRINT.ge.7.and.mod(MSTEP,IOUTPUT).eq.0) 276 + WRITE(*,'(8f10.1)')TRYCKL,TRYCKE,TRYCKB,TRYCKS,TRYCKP,TRID 277 if(IPRINT.ge.6.and.LSCTD) 278 + write(*,'(8f10.5)')(TEMPR(I),SCM(I),I=1,NTYPES) 279*** Put you own subroutine for evaluation of whatever you want here: 280*** call USER 281*** 282 end if ! if(mod(MSTEP 283* 1.4.8 Intermediate output after "NAVER" steps 284 if(mod(MSTEP,NAVER).eq.0.and.JBEAD.eq.0)CALL GETAVR(2) 285* 1.4.9 Dump restart file 286C Important - all nodes 287 if(LDMP.and.mod(MSTEP,NDUMP).eq.0)call RESTRT(2) 288 timeg = timeg+cputime(timeg0) 289* 290* 1.4.10 Conclude MD step 291* Check STOP file (visual regime only) 292 if(LVISUAL)then 293 open(unit=77,file='MD_STOP',status='old',err=10) 294 close(77) 295 return 296 end if 297 10 if(NSTEP.lt.NSTEPS)go to 1 298*<--------------------- 299* PIMD termostate analysis 300 write(*,*)' Node ',JBEAD,' averager temperature ',TEMPNOD/NSTEPS 301 if(IPRINT.ge.5)then 302 write(*,*)' particle temperatures and thermostats factors:' 303 do I=1,NSTOT 304 write(*,'(i6,f11.4,10f10.6)') 305 + I,TEMBS(I)/IPIMA,(THERMS(I,J)/(FNH*IPIMA),J=1,JNH) 306 end do 307 write(*,*)' average lengths of springs:' 308 do J=1,NSITES 309 write(*,'(a,i4,f12.4)') 310 + ' site ',J,PIRS(J)/IPIMA 311 end do 312 end if 313 RETURN 314 END 315C 316C 2. Calculate PI harminic forces between the beads 317C 318* 319*==================== PIFORCES ============================================= 320* 321 subroutine PIFORCES 322 include "prcm.h" 323 include "mpif.h" 324 real*8 PIRRT(NS) 325 VIRPI=0. 326C 327* 2.1 Retrive coordinates from upper bead (copied into array OX,,,) 328 call PI_SENDCOORD(1) 329 do I=1,NSITES 330 PIRRT(I)=0. 331 end do 332 QEK=0. 333C 334* 2.2 step forward 335 do I=1,NSTOT 336 IS = NSITE(I) 337 BX = SX(I)-OX(I) 338 BY = SY(I)-OY(I) 339 BZ = SZ(I)-OZ(I) 340 call PBC(BX,BY,BZ) 341 BSQ = BX*BX+BY*BY+BZ*BZ 342 BBB = DSQRT(BSQ) 343 R1I = 1.D0/BBB 344C Potential parameter 345 DBND = BBB*BEADFAC*MASS(IS) 346C Energy of the bond 347 EBDD = DBND*BBB 348C Force (abs. value) 349 FORB = -2.0D0*DBND*R1I 350 HX(I) = HX(I)+BX*FORB 351 HY(I) = HY(I)+BY*FORB 352 HZ(I) = HZ(I)+BZ*FORB 353 QEK = QEK + EBDD ! PI kin energy 354 PIRRT(IS)=PIRRT(IS)+BBB ! Average spring length 355 VIRPI = VIRPI + FORB*BSQ 356 VIRFX = VIRFX+FORB*BX**2*NBEADS ! will be scaled back by 1/NBEADS 357 VIRFY = VIRFY+FORB*BY**2*NBEADS 358 VIRFZ = VIRFZ+FORB*BZ**2*NBEADS 359 end do 360* 2.3. Retrive coordinates from upper bead (copied into array OX,,,) 361 call PI_SENDCOORD(-1) 362C 363 do I=1,NSTOT 364 IS = NSITE(I) 365 BX = SX(I)-OX(I) 366 BY = SY(I)-OY(I) 367 BZ = SZ(I)-OZ(I) 368 call PBC(BX,BY,BZ) 369 BSQ = BX*BX+BY*BY+BZ*BZ 370 BBB = DSQRT(BSQ) 371 R1I = 1.D0/BBB 372C Potential parameter 373 DBND = BBB*BEADFAC*MASS(IS) 374C Energy of the bond 375 EBDD = DBND*BBB 376C Force (abs. value) 377 FORB = -2.0D0*DBND*R1I 378 HX(I) = HX(I)+BX*FORB 379 HY(I) = HY(I)+BY*FORB 380 HZ(I) = HZ(I)+BZ*FORB 381 QEK = QEK + EBDD ! PI kin energy 382 end do 383 NSIT1=NSITES+1 384 PIRRT(NSIT1)=QEK 385 call MPI_ALLREDUCE(PIRRT,PIRR,NSIT1,MPI_DOUBLE_PRECISION, 386 + MPI_SUM,MPI_COMM_WORLD,ierr) 387 QEK=0.5*PIRR(NSIT1) 388 do IS=1,NSITES 389 ITYP=ITS(IS) 390 PIRR(IS)=PIRR(IS)/(NSPEC(ITYP)*NBEADS) 391 end do 392* write(*,*)' Eharm',QEK*PERMOL,JBEAD 393 return 394 end 395C 396C 3. Nose thermostate for beads 397C Note: separate thermostat for each atom !!! 398C todo: update by the chain of thermostat 399* 400*==================== BEADNOSE =============================== 401* 402 subroutine BEADNOSE 403 include "prcm.h" 404 include "mpif.h" 405 real*8 FEK(NTOT),BUF(NTOT) 406* 3.1 Compute faked kinetic energy 407 do I=1,NSTOT 408 BUF(I)=0.5*(VX(I)**2+VY(I)**2+VZ(I)**2)*MASSDI(I) 409 end do 410 call MPI_ALLREDUCE(BUF,FEK,NSTOT,MPI_DOUBLE_PRECISION, 411 + MPI_SUM,MPI_COMM_WORLD,ierr) 412 do I=1,NSTOT 413 ATEMP=FEK(I)*CONVEQ 414 DKE = ATEMP/TRTEMP-1.d0 415 THERMO(I,1)=THERMO(I,1) + DKE*DQQT*TSTEB 416 THERMS(I,1)=THERMS(I,1)+THERMO(I,1)**2 417C Nose-Hoover chain 418 if(JNH.gt.1)then 419 do J=2,JNH 420 THERMO(I,J) = THERMO(I,J) + (THERMO(I,J-1)**2-FNH)*TSTEB 421 THERMO(I,J-1) = THERMO(I,J-1)*(1.-THERMO(I,J)*TSTEB) 422 THERMS(I,J)=THERMS(I,J)+THERMO(I,J)**2 423 end do 424 end if 425 SCV=1.-THERMO(I,1)*TSTEB 426 if(IPRINT.ge.8) 427 + write(*,*)' PITH:',I,ATEMP,(THERMO(I,J),J=1,JNH),SCV 428C VX(I)=(VX(I)-FX(I))*SCV+FX(I) 429C VY(I)=(VY(I)-FY(I))*SCV+FY(I) 430C VZ(I)=(VZ(I)-FZ(I))*SCV+FZ(I) 431 VX(I)=VX(I)*SCV 432 VY(I)=VY(I)*SCV 433 VZ(I)=VZ(I)*SCV 434 TEMBS(I)=TEMBS(I)+ATEMP 435 end do 436 do I=1,NSITES 437 PIRS(I)=PIRS(I)+PIRR(I) 438 end do 439 IPIMA=IPIMA+1 440 return 441 end 442C 443C 4. Get coordinates of neighbouring node 444C ------------------------------------ 445* 446*==================== PI_SENDCOORD =============================== 447* 448 subroutine PI_SENDCOORD(IDIR) 449 include "prcm.h" 450 include "mpif.h" 451 real*4 BUFF(NBUFF),BUF2(NBUFF) 452 integer ISTS(MPI_STATUS_SIZE) 453 if(IDIR.gt.0)then 454 ISEND=JBEAD+1 455 if(ISEND.ge.NBEADS)ISEND=0 456 IRECV=JBEAD-1 457 if(IRECV.lt.0)IRECV=NBEADS-1 458 else 459 ISEND=JBEAD-1 460 if(ISEND.lt.0)ISEND=NBEADS-1 461 IRECV=JBEAD+1 462 if(IRECV.ge.NBEADS)IRECV=0 463 end if 464 NST3=3*NSTOT 465 IMA=0 466 do I=1,NSTOT 467 BUFF(I)=SX(I) 468 BUFF(I+NSTOT)=SY(I) 469 BUFF(I+2*NSTOT)=SZ(I) 470 end do 471* write(*,*)' node ',JBEAD,' in SENDCRD' 472 if(mod(JBEAD,2).eq.0)then 473 call MPI_SEND(BUFF,NST3,MPI_REAL,ISEND,IMA,MPI_COMM_WORLD,ierr) 474 call MPI_RECV(BUF2,NST3,MPI_REAL,IRECV,IMA, 475 + MPI_COMM_WORLD,ISTS,ierr) 476 else 477 call MPI_RECV(BUF2,NST3,MPI_REAL,IRECV,IMA, 478 + MPI_COMM_WORLD,ISTS,ierr) 479 call MPI_SEND(BUFF,NST3,MPI_REAL,ISEND,IMA,MPI_COMM_WORLD,ierr) 480 end if 481* call MPI_BARRIER(MPI_COMM_WORLD,IERR) 482 do I=1,NSTOT 483 OX(I)=BUF2(I) 484 OY(I)=BUF2(I+NSTOT) 485 OZ(I)=BUF2(I+2*NSTOT) 486 end do 487 return 488 end 489C 490C This compute molecular center of mass of centroids 491C 492* 493*=================== PI_COM ====================================== 494* 495 subroutine PI_COM(XPC,YPC,ZPC) 496 include "prcm.h" 497 include "mpif.h" 498 dimension XPC(NPART),YPC(NPART),ZPC(NPART) 499 double precision BUFF(NBUFF),BUF2(NBUFF) 500 do I=1,NOP 501 BUFF(I)=X(I) 502 BUFF(I+NOP)=Y(I) 503 BUFF(I+2*NOP)=Z(I) 504 end do 505 IBUF=3*NOP 506 call MPI_ALLREDUCE(BUFF,BUF2,IBUF,MPI_DOUBLE_PRECISION, 507 + MPI_SUM,MPI_COMM_WORLD,ierr) 508 do I=1,NOP 509 XPC(I)=BUF2(I)/NBEADS 510 YPC(I)=BUF2(I+NOP)/NBEADS 511 ZPC(I)=BUF2(I+2*NOP)/NBEADS 512 end do 513 return 514 end 515C 516C This compute centroid momentum by summation momentum from all nodes 517C 518* 519*=================== CMD_MOM ====================================== 520* 521 subroutine CMD_MOM 522 include "prcm.h" 523 include "mpif.h" 524 real*4 BUFF(3*NTOT),BUF2(3*NTOT) 525 do I=1,NSTOT 526 BUFF(I)=VX(I) 527 BUFF(I+NSTOT)=VY(I) 528 BUFF(I+2*NSTOT)=VZ(I) 529 end do 530 IBUF=3*NSTOT 531 call MPI_ALLREDUCE(BUFF,BUF2,IBUF,MPI_REAL, 532 + MPI_SUM,MPI_COMM_WORLD,ierr) 533 do I=1,NSTOT 534 FX(I)=BUF2(I) 535 FY(I)=BUF2(I+NSTOT) 536 FZ(I)=BUF2(I+2*NSTOT) 537 end do 538 return 539 end 540* 541*=============== SCALING ============================================== 542* 543 SUBROUTINE PI_SCALING(IAL) 544* 545* Temperature / pressure scaling 546* Path Integral version 547* ---------------------------------------- 548C 549C IAL = 1 - first half-step in double time step algorithm 550C IAL = 2 - second half-step in double time step algorithm 551C 552* 553* 1 Definitions 554* --------------- 555 include "prcm.h" 556 parameter (FP3=0.5/3.) 557 real*8 DPE(6),DPES(6),DTE(NTPS) 558 data DKEO/0./,ISKK/0/ 559* 560* 3 Recalculate temperatures and pressures 561* --------------- 562* 3.3.1 Calculate temperature 563 call GETEMP 564* 3.3.2 Calculate pressure 565C This collect virial for "atomic" pressure 566 DPE(1) = (VIR1+VIR2+VIRB+PELS1+PELS2+SPE+VIRD+ 567 + VIRA(1)+VIRA(2)+VIRA(3))/NBEADS 568C This collect virial for "molecular" pressure 569 DPE(2) = (VIR1+VIR2+PELS1+PELS2+SPE+VIRD+WIRS+WIRSS)/NBEADS 570CD write(*,'(5e13.5)')VIR1,VIR2,PELS1,PELS2,SPE,VIRD,WIRS,WIRSS 571 DPE(3) = VIRPI 572 DPE(4) = VIRFX/NBEADS 573 DPE(5) = VIRFY/NBEADS 574 DPE(6) = VIRFZ/NBEADS 575 call ALLSUM(DPE,DPES,6) 576 VIRPI = DPES(3) 577 VIRSUM = DPES(1)+VIRPI 578 WIRSUM = DPES(2)+VIRPI 579 VIRFX = DPES(4) 580 VIRFY = DPES(5) 581 VIRFZ = DPES(6) 582C Kinetic (ideal gas) contributions to pressure 583 TRID = 2.*TKE*UNITP*NBEADS/(3.*VOL) 584 TRIDM = NOP*BOLTZ*TEMP*1.d-5*NBEADS/(VOL*UNITL**3) 585C Pressure in atm and its projections 586 TRYCK = VIRSUM*UNITP/(3.*VOL)+TRID 587* TRYCKM = WIRSUM*UNITP/(3.*VOL)+TRIDM 588 TRYCKM=0. 589C Projections are calculated for "atomic" pressure 590 TRYCKX = (VIRX+VIRFX+VIRD/3.+VIRA(1))*UNITP/VOL+TRID 591 TRYCKY = (VIRY+VIRFY+VIRD/3.+VIRA(2))*UNITP/VOL+TRID 592 TRYCKZ = (VIRZ+VIRFZ+VIRD/3.+VIRA(3))*UNITP/VOL+TRID 593CD if(TASKID.eq.MAST)write(*,'(5f13.2)') 594CD +TRYCK,TRYCKX,TRYCKY,TRYCKZ,(TRYCKX+TRYCKY+TRYCKZ)/3. 595C 596C Now temperature TEMP and pressure TRYCK are defined for the 597C current configuration, as well as their components 598* 599* 3.3.3 Calculate deviation from thermostat T and P 600 DKE = TEMP/TRTEMP-1.d0 ! total temperature 601 do I=1,NTYPES 602CD write(*,*)' temp ',I,TEMPR(I),TASKID 603 DTE(I) = TEMPR(I)/TRTEMP-1.d0 ! for each species 604 end do 605C in the case of constrained dynamics, pressure is defined by 606C "molecule" algorithm - it corresponds scaling of molecular COM 607C for flexible molecules, scaling of atom positions is employed, 608C which corresponds to the "atomic" pressure 609C Exception: case of separate pressure control in each direction 610C (LSEP=.true.). Then pressure in each direction is determined from 611C "atomic" pressure. 612 if(LSHEJK)then 613 DPE(1) = TRYCKM-TRPRES 614 else 615 DPE(1) = TRYCK -TRPRES 616 end if 617 DPE(2) = TRYCKX-TRPRES 618 DPE(3) = TRYCKY-TRPRES 619 if(LHEX)then 620 DPE(2)=0.5*(DPE(2)+DPE(3)) 621 DPE(3)=DPE(2) 622 end if 623 DPE(4) = TRYCKZ-TRPRES 624* 625* 3.4 First integration of Nose-Hoover equations (t->t+1/2) 626* ----------------------------------------------------------- 627 if(IAL.eq.1)then 628 if(LNVT)then 629* 3.4.1 Correction ksi due to thermostat 630C coefficients DQT,DQP and RTP were defined in main.f, p.6.2 631C 632C Nose eqn is: dP/dt -> F - P*ksi/Q 633C dksi/dt -> 2*Ekin - Nf*kT 634C Here ksi/Q = SC 635C Q = Nf*kT*tau**2 ; tau -> QT (i.u.) 636C Q/(Nf*kT) -> 1./DQT 637C So d(SC)/DT = (2*Ekin/Nf*kT-1)*DQT = DKE*DQT 638C 639C This is separate thermostat for each species 640 if(LSCTD)then 641 SCA = 0. 642 do I=1,NTYPES 643 SCM(I) = SCM(I) + DTE(I)*DQT*HTSTEP ! DQT -thermostat mass 644 SCA = SCA + SCM(I)*EKIN(I) 645 end do 646 SC = SCA/TKE 647C common thermostat 648 else 649 SC = SC + DKE*DQT*HTSTEP 650 end if 651 if(LNPT)then 652* 3.4.2 Barostat corrections (NPT) 653* ------------------------------- 654* 3.4,2,1 Correction ksi due to barostat 655C 656C Additional correction to ksi: 657C d(ksi)=ita**2/W-kT 658C ita/W -> SCL 659C so d(ksi) -> SCL**2*(W/Q) - kT/Q 660C 661 if(LSCTD)then 662 SCA = 0. 663 do I=1,NTYPES 664 SCM(I) = SCM(I) + (SCL**2*RTP-DQT/FNST)*HTSTEP 665 SCA = SCA + SCM(I)*EKIN(I) 666 end do 667 SC = SCA/TKE 668 else 669 SC = SC + (SCL**2*RTP-DQT/FNST)*HTSTEP ! RTP=DQT/DQP 670 end if 671* 3.4.2.2 Correction ita 672C TKE is kinetic energy 673 if(LSEP)then 674 SCLX = SCLX+DQP*(VOL*DPE(2)+2.*TKE/FNST-SCLX*SC)*HTSTEP 675 SCLY = SCLY+DQP*(VOL*DPE(3)+2.*TKE/FNST-SCLY*SC)*HTSTEP 676 if(LHEX)then 677 SCLX=0.5*(SCLX+SCLY) 678 SCLY=SCLX 679 end if 680 SCLZ = SCLZ+DQP*(VOL*DPE(4)+2.*TKE/FNST-SCLZ*SC)*HTSTEP 681 SCL = SCLX+SCLY+SCLZ ! trace(P) 682 else 683 SCL = SCL + DQP*(3.*VOL*DPE(1)+6.*TKE/FNST-SCL*SC)*HTSTEP 684 SCLX=SCL 685 SCLY=SCL 686 SCLZ=SCL 687C Control fluctuations 688 if(abs(SCL*HTSTEP).ge.0.1)then 689 write(*,*)' too strong fluctuations in NPT algorithm' 690 if(NUMTASK.gt.1)write(*,*)' at node ',TASKID 691 write(*,*)' scaling factor ',SCL*HTSTEP 692 if(IPRINT.ge.7)then 693 write(*,*)'------------------->' 694 write(*,*)SCL,TRYCK,TRYCKM 695 write(*,*)VIR1*UNITP/(3.*VOL),VIR2*UNITP/(3.*VOL), 696 + PELS1*UNITP/(3.*VOL),PELS2*UNITP/(3.*VOL),VIRB*UNITP/(3.*VOL) 697 write(*,*)'<-------------------' 698 end if 699 ISKK=ISKK+1 700 if(ISKK.gt.100)then 701 write(*,*)' !!! repeated failure of NPT-algorithm' 702 write(*,*)' !!! restart program with constant volume' 703 write(*,*)' !!! or increase thermostat parameter for pressure' 704 call FINAL 705 end if 706 end if ! if(LSEP 707 end if 708* 3.4.2.3 Calculate scaling coefficients 709 SCVX = 1.-((1.-3./FNST)*SCL+SC)*HTSTEP 710 SCVY = SCVX 711 SCVZ = SCVX 712 DRCX = TSTEP*SCLX 713 DRCY = TSTEP*SCLY 714 DRCZ = TSTEP*SCLZ 715 DRC = (DRCX+DRCY+DRCZ)/3. 716 DRCV = DRC*HTSTEP ! These are second order corrections 717 DRCX = 1.+DRCX + 0.5*DRCX**2 ! and (may be?) ommitted 718 DRCY = 1.+DRCY + 0.5*DRCY**2 ! and (may be?) ommitted 719 DRCZ = 1.+DRCZ + 0.5*DRCZ**2 ! and (may be?) ommitted 720* 721* 3.4.2.4 Correct coordinates and velocities - flexible molecules 722 do I=NAB(TASKID),NAE(TASKID) 723 ITYP = ITYPE(I) 724 if(LMOVE(ITYP))then 725 VX(I) = VX(I)*SCVX 726 VY(I) = VY(I)*SCVY 727 VZ(I) = VZ(I)*SCVZ 728 SX(I) = SX(I)*DRCX + VX(I)*DRCV*MASSDI(I) 729 SY(I) = SY(I)*DRCY + VY(I)*DRCV*MASSDI(I) 730 SZ(I) = SZ(I)*DRCZ + VZ(I)*DRCV*MASSDI(I) 731 end if 732 end do 733* 3.4.2.5 Scale simulation box 734 call RECLEN(DRCX,DRCY,DRCZ) 735 else ! .not.LNPT 736* 3.4.3 Velocities corrections in the NVT ensemble 737 call GETEMP 738 if(LSCTD)then 739 do I=NAB(TASKID),NAE(TASKID) 740 ITYP = ITYPE(I) 741 SCV = 1.-SCM(ITYP)*HTSTEP 742 if(LMOVE(ITYP))then 743 VX(I) = VX(I)*SCV 744 VY(I) = VY(I)*SCV 745 VZ(I) = VZ(I)*SCV 746 end if 747 end do 748 else 749 SCV = 1.-SC*HTSTEP 750 do I=NAB(TASKID),NAE(TASKID) 751 ITYP = ITYPE(I) 752 if(LMOVE(ITYP))then 753 VX(I) = VX(I)*SCV 754 VY(I) = VY(I)*SCV 755 VZ(I) = VZ(I)*SCV 756 end if 757 end do 758 end if ! if(LSCTD) 759 end if ! if(LNPT) 760* 3.4.4 Absorbed energy 761 EABS = EABS - (SC+SCL*(1.-3./FNST))*TKE*TSTEP 762 end if ! if(LNVT) 763* 764* 3.5 Second integration of Nose-Hoover equations (t+1/2->t+1) 765* ------------------------------------------------------------- 766 else if(IAL.eq.2)then 767 if(LNVT)then 768* 3.5.1 Correct velocities due to thermostat 769 if(LSCTD)then 770 do I=NAB(TASKID),NAE(TASKID) 771 ITYP = ITYPE(I) 772 SCV = 1.-SCM(ITYP)*HTSTEP 773 if(LMOVE(ITYP))then 774 VX(I) = VX(I)*SCV 775 VY(I) = VY(I)*SCV 776 VZ(I) = VZ(I)*SCV 777 end if 778 end do 779 else 780 SCV = 1.-SC*HTSTEP 781 do I=NAB(TASKID),NAE(TASKID) 782 ITYP = ITYPE(I) 783 if(LMOVE(ITYP))then 784 VX(I) = VX(I)*SCV 785 VY(I) = VY(I)*SCV 786 VZ(I) = VZ(I)*SCV 787 end if 788 end do 789 end if 790* 3.5.2 Correct velocities due to barostat 791 if(LNPT)then 792 SCVX = 1.-(1.-3./FNOP)*SCL*HTSTEP 793 SCVY = SCVX 794 SCVZ = SCVX 795 do I=NAB(TASKID),NAE(TASKID) 796 ITYP = ITYPE(I) 797 if(LMOVE(ITYP))then 798 VX(I) = VX(I)*SCVX 799 VY(I) = VY(I)*SCVY 800 VZ(I) = VZ(I)*SCVZ 801 end if 802 end do 803* 3.5.3 Recalculate temperature 804 call GETEMP 805* 3.5.4 Correct ita 806 if(LSEP)then 807 SCLX = SCLX+DQP*(VOL*DPE(2)+2.*TKE/FNST-SCLX*SC)*HTSTEP 808 SCLY = SCLY+DQP*(VOL*DPE(3)+2.*TKE/FNST-SCLY*SC)*HTSTEP 809 if(LHEX)then 810 SCLX=0.5*(SCLX+SCLY) 811 SCLY=SCLX 812 end if 813 SCLZ = SCLZ+DQP*(VOL*DPE(4)+2.*TKE/FNST-SCLZ*SC)*HTSTEP 814 SCL = SCLX+SCLY+SCLZ ! trace(P) 815 else 816 SCL = SCL + DQP*(3.*VOL*DPE(1)+6.*TKE/FNST-SCL*SC)*HTSTEP 817 SCLX=SCL 818 SCLY=SCL 819 SCLZ=SCL 820C Control fluctuations 821 if(abs(SCL*HTSTEP).ge.0.1)then 822 write(*,*)' too strong fluctuations in NPT algorithm' 823 if(NUMTASK.gt.1)write(*,*)' at node ',TASKID 824 write(*,*)' scaling factor ',SCL*HTSTEP 825 if(IPRINT.ge.7)then 826 write(*,*)'------------------->' 827 write(*,*)SCL,TRYCK,TRYCKM 828 write(*,*)VIR1*UNITP/(3.*VOL),VIR2*UNITP/(3.*VOL), 829 + PELS1*UNITP/(3.*VOL),PELS2*UNITP/(3.*VOL),VIRB*UNITP/(3.*VOL) 830 write(*,*)'<-------------------' 831 end if 832 ISKK=ISKK+1 833 if(ISKK.gt.100)then 834 write(*,*)' !!! repeated failure of NPT-algorithm' 835 write(*,*)' !!! restart program with constant volume' 836 write(*,*)' !!! or increase thermostat parameter for pressure' 837 call FINAL 838 end if 839 end if ! if(LSEP 840 end if 841* 3.5.5 Correct ksi due to barostat 842 SC = SC + (SCL**2*RTP-DQT/FNST)*HTSTEP 843 if(LSCTD)then 844 do I=1,NTYPES 845 SCM(I) = SCM(I) + (SCL**2*RTP-DQT/FNST)*HTSTEP 846 end do 847 end if 848 else ! .not.LNPT 849* recalculate temperature 850 call GETEMP 851 end if ! if(LNPT 852* 3.5.6 Correction ksi due to thermostat 853 if(LSCTD)then 854 SCA = 0. 855 do I=1,NTYPES 856 DTE(I) = TEMPR(I)/TRTEMP-1.d0 857 SCM(I) = SCM(I) + DTE(I)*DQT*HTSTEP 858 SCA = SCA + SCM(I)*EKIN(I) 859 end do 860 SC = SCA/TKE 861 else 862 DKE = TEMP/TRTEMP-1.d0 863 SC = SC + DKE*DQT*HTSTEP 864 end if 865* 3.5.7 Absorbed energy 866 EABS = EABS - (SC+SCL*(1.-3./FNST))*TKE*TSTEP 867 end if ! if (LNVT 868* 869* 3.6 Case of single-step SHAKE algorithm 870* ---------------------------------------- 871 else if(IAL.eq.3)then 872* 3.6.1 Correct ksi(t+1) from velocities and temperature at t+1/2 873 if(LNVT)then 874 SC = SC + DKE*DQT*TSTEP ! Ksi(t+1) 875C This is separate thermostat for each species (may be not work) 876 if(LSCTD)then 877 SCA = 0. 878 do I=1,NTYPES 879 SCM(I) = SCM(I) + DTE(I)*DQT*TSTEP ! DQT -thermostat mass 880 SCA = SCA + SCM(I)*EKIN(I) 881 end do 882 SC = SCA/TKE 883 end if 884 if(LNPT)then 885* 3.6.2 Barostat corrections (NPT) 886* ------------------------------- 887* 3.6.2.1 Correction ita t-1/2 -> t+1/2 888C TKE is kinetic energy 889 if(LSEP)then 890 SCLX = SCLX+DQP*(VOL*DPE(2)+2.*TKE/FNST-SCLX*SC)*TSTEP 891 SCLY = SCLY+DQP*(VOL*DPE(3)+2.*TKE/FNST-SCLY*SC)*TSTEP 892 if(LHEX)then 893 SCLX=0.5*(SCLX+SCLY) 894 SCLY=SCLX 895 end if 896 SCLZ = SCLZ+DQP*(VOL*DPE(4)+2.*TKE/FNST-SCLZ*SC)*TSTEP 897 SCL = SCLX+SCLY+SCLZ ! trace(P) 898 else 899 SCL = SCL + DQP*(3.*VOL*DPE(1)+6.*TKE/FNST-SCL*SC)*TSTEP 900 SCLX=SCL 901 SCLY=SCL 902 SCLZ=SCL 903* 3.6.2.2 Control fluctuations 904 if(abs(SCL*TSTEP).ge.0.1)then 905 write(*,*)' too strong fluctuations in NPT algorithm' 906 if(NUMTASK.gt.1)write(*,*)' at node ',TASKID 907 write(*,*)' scaling factor ',SCL*TSTEP 908 if(IPRINT.ge.7)then 909 write(*,*)'------------------->' 910 write(*,*)SCL,TRYCK,TRYCKM 911 write(*,*)VIR1*UNITP/(3.*VOL),VIR2*UNITP/(3.*VOL), 912 + PELS1*UNITP/(3.*VOL),PELS2*UNITP/(3.*VOL),VIRB*UNITP/(3.*VOL) 913 write(*,*)'<-------------------' 914 end if 915 ISKK=ISKK+1 916 if(ISKK.gt.100)then 917 write(*,*)' !!! repeated failure of NPT-algorithm' 918 write(*,*)' !!! restart program with constant volume' 919 write(*,*)' !!! or increase thermostat parameter for pressure' 920 call FINAL 921 end if 922 end if ! if(LSEP 923 end if 924* this is for correction of velocities 925 SCV = ((1.-3./FNST)*SCL)*TSTEP 926* 3.6.2.3 Calculate scaling coefficients 927 DRCX = TSTEP*SCLX ! eta*dt 928 DRCY = TSTEP*SCLY 929 DRCZ = TSTEP*SCLZ 930* 3.6.2.4 Correct coordinates and velosities 931* 3.6.2.4.1 Calculate centre of mass and molecular momenta 932 IMB = NNUM(NAB(TASKID)) 933 IME = NNUM(NAE(TASKID)) 934C Cycle over molecules 935 do IMOL = IMB,IME 936 ITYP = ITM(IMOL) 937 SUMM = SUMMAS (ITYP) 938 if(LMOVE(ITYP))then 939 NSBEG = ISADR(ITYP)+1 940 NSEND = ISADR(ITYP+1) 941 XC = 0.D0 942 YC = 0.D0 943 ZC = 0.D0 944 PXC = 0. 945 PYC = 0. 946 PZC = 0. 947C Cycle over atoms within the molecule 948 IBEG=ISADDR(ITYP)+1+NSITS(ITYP)*(IMOL-IADDR(ITYP)-1) 949 IEND=ISADDR(ITYP)+NSITS(ITYP)*(IMOL-IADDR(ITYP)) 950 DO I = IBEG,IEND 951 IS = NSITE(I) 952 XC = XC+MASS(IS)*SX(I) 953 YC = YC+MASS(IS)*SY(I) 954 ZC = ZC+MASS(IS)*SZ(I) 955C Molecular momenta 956 PXC = PXC+VX(I) 957 PYC = PYC+VY(I) 958 PZC = PZC+VZ(I) 959 END DO! OF N 960 XC = XC/SUMM 961 YC = YC/SUMM 962 ZC = ZC/SUMM 963* 3.6.2.4.2 Correction to COM momenta 964 DVX = PXC*SCV 965 DVY = PYC*SCV 966 DVZ = PZC*SCV 967* 3.6.2.4.3 COM displacements 968 DXC = XC*DRCX 969 DYC = YC*DRCY 970 DZC = ZC*DRCZ 971* 3.6.2.4.4 Correct atom velocities 972C DVX - molecular momentum correction 973C DVX/SUMM - COM velocity correction 974C MASS(IS)*DVX/SUMM - atom momentum correction 975 DO I = IBEG,IEND 976 IS = NSITE(I) 977 FAC = MASS(IS)/SUMM 978 VX(I) = VX(I)-DVX*FAC 979 VY(I) = VY(I)-DVY*FAC 980 VZ(I) = VZ(I)-DVZ*FAC 981* 3.6.2.4.5 Correct atom positions 982 SX(I) = SX(I) + DXC 983 SY(I) = SY(I) + DYC 984 SZ(I) = SZ(I) + DZC 985 end do 986 end if 987 end do ! of ITYP 988* 3.6.2.5 Scale simulation box 989 DRCX=DRCX+1. 990 DRCY=DRCY+1. 991 DRCZ=DRCZ+1. 992 call RECLEN(DRCX,DRCY,DRCZ) 993* 3,6,2,6 Correction ksi due to barostat 994C 995C Additional correction to ksi: 996C d(ksi)=ita**2/W-kT 997C ita/W -> SCL 998C so d(ksi) -> SCL**2*(W/Q) - kT/Q 999C 1000 SC = SC + (SCL**2*RTP-DQT/FNST)*TSTEP ! RTP=DQT/DQP 1001 if(LSCTD)then 1002 SCA = 0. 1003 do I=1,NTYPES 1004 SCM(I) = SCM(I) + (SCL**2*RTP-DQT/FNST)*TSTEP 1005 SCA = SCA + SCM(I)*EKIN(I) 1006 end do 1007 SC = SCA/TKE 1008 end if 1009 end if ! if(LNPT) 1010* 3.6.4 Absorbed energy 1011 EABS = EABS - (SC+SCL*(1.-3./FNST))*TKE*TSTEP 1012 end if ! if(LNVT) 1013 end if ! if(IAL 1014* 1015* 3.7 Forcible adjusting of velocities to given temperature 1016* --------------------------------------------------------- 1017 IF(LSCLT.and.(IAL.eq.2.or.IAL.eq.3)) THEN 1018 if(LSCTD)then 1019 DO ITYP = 1,NTYPES 1020 IF(DABS(TRTEMP-TEMPR(ITYP)).GT.TDELT.and.LMOVE(ITYP))THEN 1021 SCV = DSQRT(TRTEMP/TEMPR(ITYP)) 1022 IBEG=ISADDR(ITYP)+1 1023 IEND=ISADDR(ITYP+1) 1024 if(IBEG.le.NAB(TASKID))IBEG=NAB(TASKID) 1025 if(IEND.gt.NAE(TASKID))IEND=NAE(TASKID) 1026 do I=IBEG,IEND 1027 VX(I) = SCV*VX(I) 1028 VY(I) = SCV*VY(I) 1029 VZ(I) = SCV*VZ(I) 1030 end do 1031 if(IPRINT.ge.5)write(*,*) 1032 + 'velocities scaled by ',SCV,' for type',ITYP 1033 NRTSC(ITYP) = NRTSC(ITYP)+1 1034 end if 1035 end do 1036 else 1037 IF(DABS(TRTEMP-TEMP).GT.TDELT)THEN 1038 SCV = dsqrt(TRTEMP/TEMP) 1039 do I=NAB(TASKID),NAE(TASKID) 1040 ITYP = ITYPE(I) 1041 if(LMOVE(ITYP))then 1042 VX(I) = VX(I)*SCV 1043 VY(I) = VY(I)*SCV 1044 VZ(I) = VZ(I)*SCV 1045 end if 1046 end do 1047 if(IPRINT.ge.5)write(*,*) 1048 + 'velocities scaled by ',SCV 1049 NRTSC(ITYP) = NRTSC(ITYP)+1 1050 end if 1051 end if 1052 call GETEMP 1053 END IF!(LSC... 1054* 1055 RETURN 1056 END 1057* 1058*================== PI_AVER ====================================== 1059* 1060 subroutine PI_AVER 1061 include "prcm.h" 1062 include "mpif.h" 1063 real*8 BUFF(NBUFF),BUF2(NBUFF) 1064 BUFF(1) = PELS1 1065 BUFF(2) = PELS2 1066 BUFF(3) = SPE 1067 BUFF(4) = PE1 1068 BUFF(5) = PE2 1069 BUFF(6) = PEST 1070 BUFF(7) = PINT 1071 BUFF(8) = TKE 1072 BUFF(9) = TROT 1073 BUFF(10) = TTR 1074 BUFF(11) = TINT 1075 BUFF(12) = TEMP 1076 BUFF(13)=VIR1 1077 BUFF(14)=VIR2 1078 BUFF(15)=VIRB 1079 BUFF(16)=WIRSUM 1080 BUFF(17) = VIRD 1081 IBUF=17+MOLINT 1082 do INT=1,MOLINT 1083 BUFF(17+INT)=POTLJ(INT) 1084 BUFF(IBUF+INT)=POTES(INT) 1085 end do 1086 IBUF=IBUF+MOLINT 1087 IBUF2=IBUF+NTYPES 1088 do I=1,NTYPES 1089 BUFF(IBUF+I)=PES14(I) 1090 BUFF(IBUF2+I)=PSR14(I) 1091 end do 1092 IBUF=IBUF2+NTYPES 1093 if(ISTAV.gt.0)then 1094 if(NRBB.gt.0)then 1095 do I=1,NRBB 1096 BUFF(IBUF+I)=BE(I) 1097 end do 1098 IBUF=IBUF+NRBB 1099 end if 1100 if(NRAA.gt.0)then 1101 do I=1,NRAA 1102 BUFF(IBUF+I)=AE(I) 1103 end do 1104 IBUF=IBUF+NRAA 1105 end if 1106 if(NRTT.gt.0)then 1107 do I=1,NRTT 1108 BUFF(IBUF+I)=TE(I) 1109 end do 1110 IBUF=IBUF+NRTT 1111 end if 1112 end if ! if(ISTAV 1113 if(IBUF.gt.NBUFF)then 1114 write(*,*) 'increase NBUFF to ',IBUF 1115 call FINAL 1116 end if 1117 call MPI_ALLREDUCE(BUFF,BUF2,IBUF,MPI_DOUBLE_PRECISION, 1118 + MPI_SUM,MPI_COMM_WORLD,ierr) 1119 PELS1 = BUF2(1)/NBEADS 1120 PELS2 = BUF2(2)/NBEADS 1121 SPE = BUF2(3)/NBEADS 1122 PE1 = BUF2(4)/NBEADS 1123 PE2 = BUF2(5)/NBEADS 1124 PEST = BUF2(6)/NBEADS 1125 PINT = BUF2(7)/NBEADS 1126 TKE = BUF2(8)/NBEADS 1127 TROR = BUF2(9)/NBEADS 1128 TTR = BUF2(10)/NBEADS 1129 TINT = BUF2(11)/NBEADS 1130 TEMP = BUF2(12)/NBEADS 1131 VIR1 = BUF2(13)/NBEADS 1132 VIR2 = BUF2(14)/NBEADS 1133 VIRB = BUF2(15)/NBEADS 1134 WIRSUM = BUF2(16)/NBEADS 1135 VIRD = BUF2(17)/NBEADS 1136 IBUF=17+MOLINT 1137 do INT=1,MOLINT 1138 POTLJ(INT)=BUF2(INT+17)/NBEADS 1139 POTES(INT)=BUF2(IBUF+INT)/NBEADS 1140 end do 1141 IBUF=IBUF+MOLINT 1142 IBUF2=IBUF+NTYPES 1143 do I=1,NTYPES 1144 PES14(I)=BUF2(IBUF+I)/NBEADS 1145 PSR14(I)=BUF2(IBUF2+I)/NBEADS 1146 end do 1147 IBUF=IBUF2+NTYPES 1148 if(ISTAV.gt.0)then 1149 if(NRBB.gt.0)then 1150 do I=1,NRBB 1151 BE(I)=BUF2(IBUF+I)/NBEADS 1152 end do 1153 IBUF=IBUF+NRBB 1154 end if 1155 if(NRAA.gt.0)then 1156 do I=1,NRAA 1157 AE(I)=BUF2(IBUF+I)/NBEADS 1158 end do 1159 IBUF=IBUF+NRAA 1160 end if 1161 if(NRTT.gt.0)then 1162 do I=1,NRTT 1163 TE(I)=BUF2(IBUF+I)/NBEADS 1164 end do 1165 IBUF=IBUF+NRTT 1166 end if 1167 end if 1168 return 1169 end 1170