1* Program MDynaMix, v.5.2 2* 3* Part 5 4* 5* File forces.f 6* ------------- 7* 8C This file contains subroutines responsible for calculations 9C of different forces. Forces are divided into two groups: 10C slow and fast forces. This division is essential only if 11C double time step algorithm is used 12C 13C Subroutines: 14C 15C 1. FFORCES - collect fast forces 16C 2. SLFORCES - collect slow forces (double time step) 17C 3. ALLFORCE - collect all forces (used in simple verlet) 18C 19C 4. GETBND - calculate covalent bonds forces 20C 5. GETANG - calculate covalent angle forces 21C 6. GETTRS - calculate forces due to dihedral angles (all types) 22C 7. STRBND - a special procedure for SPC water model 23C 24C 8. LOCALF - non-bonded forces of closest neighbours 25C 9. FORCES - non-bonded forces of far neighbours 26C 10. FURIR - reciprical space Ewald contribution 27C 11. INIFUR - initialisation of FURIR subroutine 28C 12. ETERMS - self-interaction term of Ewald sum and intramolecular corrections 29C 13. ELRCLJ - long-range corrections from LJ forces 30C 14. PULLBACK - fixed harmonic potential for selected atoms 31 32C 15. CHCNB - update list of neighbours 33* 34* 1. Collect fast forces 35* ---------------------- 36 subroutine FFORCES(LAVER) 37 logical LAVER 38* 1.1 Set zero-s 39 CALL ZEROFS(1) ! service.f 40* 1.2 Calculate different contributions to the fast forces 41C (see details in descriptions of corresponding subroutines below) 42 call LOCALF 43 CALL GETBND 44 CALL GETANG 45 CALL GETTRS 46* 1.3 Sum up fast forces from different nodes 47 call CM_ADDF(LAVER) ! mpi.f or scalar.f 48* 1.4 Cut large forces of specified 49 call SCLFRC(1) 50 return 51 end 52* 53* 2. Collect slow forces 54* ---------------------- 55 subroutine SLFORCES 56* 2.1 Set zero-s 57 call ZEROAV ! service.f 58 CALL ZEROFS(2) ! service.f 59* 2.2 Calculate different contributions to the slow forces 60C (see details in descriptions of corresponding subroutines below) 61 CALL FURIR ! 62 CALL ETERMS ! 63 CALL FORCES ! l-forces.f 64* 2.3 Sum up "slow" forces from different nodes 65 CALL CM_ADDLF ! mpi.f 66C This may take essential time but seems has no effect 67* if(LNPT)call ELRCLJ ! l-forces.f 68C If a molecule is softly fixed 69 call PULLBACK 70* 2.4 Cut large forces of specified 71 call SCLFRC(2) 72 return 73 end 74* 75* 3. Collect all forces 76* ---------------------- 77 subroutine ALLFORCE 78* 3.1 Set zero-s 79 call ZEROAV ! service.f 80 CALL ZEROFS(1) ! service.f 81 CALL ZEROFS(2) ! service.f 82* 3.2 Calculate different contributions to the forces 83C (see details in descriptions of corresponding subroutines below) 84 call LOCALF 85 CALL GETBND 86 CALL GETANG 87 CALL GETTRS 88 CALL FURIR 89 CALL ETERMS 90 CALL FORCES 91* 3.3 Sum up all forces from different nodes 92 CALL CM_ADDLA ! mpi.f 93C This may take essential time but seems has no effect 94* if(LNPT)call ELRCLJ ! l-forces.f 95C If a molecule is softly fixed 96 call PULLBACK 97* 3.4 Cut large forces of specified 98 call SCLFRC(3) 99 return 100 end 101* 102*=============== INTRAMOLECUULAR FORCES ========================== 103* 104* 4. Covalent bond contribution 105* ----------------------------- 106C 107C This subroutine calculate contribution to forces from 108C covalent bonds. Each node has own list of bonds. Distribution 109C of bonds overthe nodes takes place in subroutine UNITS (setup.f) 110C 111C=============== GETBND =========================================== 112C 113 SUBROUTINE GETBND 114* 115* 4.1 Front matter 116* ---------------- 117 include "prcm.h" 118 timeb0 = cputime(0.d0) 119C This set zeroth in arrays calculation average bond 120 DO M = 1,NRBB 121 BR(M) = 0.D0 122 BE(M) = 0.D0 123 END DO 124C This call special procedure to calculate bonds in the case if 125C parameter IPOT for a molecule is 1 or 2 126 do ITYP=1,NTYPES 127 if(IPOT(ITYP).eq.1.or.IPOT(ITYP).eq.2)call STRBND(ITYP) 128 end do 129 if(NRBB.le.0)go to 200 130* 131* 4.2 Organize cycle over bonds and define parameters of each bond 132* ---------------------------------------------------------------- 133C Cycle over the bonds for this node 134 do N=1,NRBT 135 I = IBI(N) ! first atom 136 ITYP = ITYPE(I) ! type 137 M = IBUK(N) ! global bond number 138 J = JBJ(N) ! second atom 139 NSP = NSPEC(ITYP) ! num. of molecules of this type 140 FNSPI = 1.D0/DFLOAT(NSP) 141C Equilibrium length of the bond 142 REQ = RB(M) 143* 144* 4.3 Calculate bond length 145* ------------------------- 146 BX = SX(J)-SX(I) 147 BY = SY(J)-SY(I) 148 BZ = SZ(J)-SZ(I) 149 call PBC(BX,BY,BZ) 150 BSQ = BX*BX+BY*BY+BZ*BZ 151 BD = DSQRT(BSQ) 152 R1I = 1.D0/BD 153C Deviation fron equilibrium 154 BBB = BD-REQ 155* 156* 4.4 Calculate energy and force 157* ------------------------------ 158* 4.4.1 Case of harmonic bond 159 if(ID(M).ne.1.or.BBB.gt.REQ)then 160C Potential parameter 161 FCN = FB(M) 162 DBND = BBB*FCN 163C Energy of the bond 164 EBDD = DBND*BBB 165C Force (abs. value) 166 FORB = -2.0D0*DBND*R1I 167* 4.4.2 Case of Morse potential 168 else 169 DDIS = DB(M) 170 ROH = RO(M) 171 EXM = exp(-ROH*BBB) 172 EXMP = 1-EXM 173C Energy of the bond 174 EBDD = DDIS*EXMP**2 175C Force (abs. value) 176 FORB = -2.d0*DDIS*ROH*EXM*EXMP*R1I 177 end if 178* 179* 4.5 Add results to accumullating arrays 180* --------------------------------------- 181 BR(M) = BR(M)+BD*FNSPI ! average length 182 BE(M) = BE(M)+EBDD*FNSPI ! average energy 183 PINT = PINT+EBDD ! total intramolec. energy 184 if(LMOVE(ITYP))then 185C Forces 186 HX(I) = HX(I)-BX*FORB 187 HY(I) = HY(I)-BY*FORB 188 HZ(I) = HZ(I)-BZ*FORB 189 HX(J) = HX(J)+BX*FORB 190 HY(J) = HY(J)+BY*FORB 191 HZ(J) = HZ(J)+BZ*FORB 192C Contribution to the virial 193 VIRB = VIRB+FORB*BSQ 194 VIRFX = VIRFX+FORB*BX**2 195 VIRFY = VIRFY+FORB*BY**2 196 VIRFZ = VIRFZ+FORB*BZ**2 197 end if 198 if(IPRINT.ge.9)write(*,*)I,J,BD,EBDD/EFACT 199* 200 END DO! OF M 201* 202 200 timeb = timeb+cputime(timeb0) 203 RETURN 204 END 205* 206* 5. Covalent angle contribution 207* ------------------------------- 208C 209C This subroutine calculate contribution to forces from 210C covalent angles. Each node has own list of angles. Distribution 211C of angles over the nodes takes place in subroutine UNITS (setup.f) 212C 213C 214C=============== GETANG =========================================== 215C 216 SUBROUTINE GETANG 217* 218* 5.1 Front matter 219* ---------------- 220 include "prcm.h" 221 IF(NRAA.LE.0) RETURN 222 timea0 = cputime(0.d0) 223* 224 DO M = 1,NRAA 225 AR(M) = 0.D0 226 AE(M) = 0.D0 227 END DO 228* 229* 5.2 Organize cycle over the angles 230* ---------------------------------- 231 do N=1,NRAT 232 I=IAI(N) ! first atom 233 ITYP=ITYPE(I) ! type 234 J = JAJ(N) ! second atom 235 K = KAK(N) ! third atom 236 M = IAUK(N) ! global angle number 237 FCN = FA(M) ! force constant 238 AEQ = RA(M) ! equilibrium angle 239 NSP = NSPEC (ITYP) 240 FNSPI = 1.D0/DFLOAT(NSP) 241* 242* 5.3 Calculate the angle and energy 243* ---------------------------------- 244* 5.3.1 I->J vector 245 AX = SX(I)-SX(J) 246 AY = SY(I)-SY(J) 247 AZ = SZ(I)-SZ(J) 248 call PBC(AX,AY,AZ) 249 ASQ = AX*AX+AY*AY+AZ*AZ 250 ASQI = 1.D0/ASQ 251* 5.3.2 K->J vector 252 BX = SX(K)-SX(J) 253 BY = SY(K)-SY(J) 254 BZ = SZ(K)-SZ(J) 255 call PBC(BX,BY,BZ) 256 BSQ = BX*BX+BY*BY+BZ*BZ 257 BSQI = 1.D0/BSQ 258* 5.3.3 Calculate Angle 259 AB = DSQRT(ASQ*BSQ) 260 ABI = 1.D0/AB 261 COSA = AX*BX+AY*BY+AZ*BZ 262 COSA = COSA*ABI 263 COSA = DMIN1(COSA, 1.D0) 264 COSA = DMAX1(COSA,-1.D0) 265 ALF = DACOS(COSA) 266 SINA = DSIN(ALF) 267 SINAI = 1.D0/SINA 268* 5.3.4 Calculate energy 269 DALF = ALF-AEQ ! deviation from equilibrium 270 ABC = DALF*FCN ! derivative of the energy 271 EANG = ABC*DALF ! energy 272 AR(M) = AR(M)+ALF*FNSPI 273 AE(M) = AE(M)+EANG*FNSPI 274 PINT = PINT+EANG 275 if(IPRINT.ge.9)write(*,*)I,J,K,ALF*TODGR,EANG/EFACT 276* 277* 5.4 Force calculation 278* --------------------- 279 if(LMOVE(ITYP))then 280 SAB =-2.D0*ABC*SINAI 281 CAB = SAB*COSA 282 FAB = SAB*ABI 283 FAA = CAB*ASQI 284 FBB = CAB*BSQI 285 FAX = FAB*BX-FAA*AX 286 FAY = FAB*BY-FAA*AY 287 FAZ = FAB*BZ-FAA*AZ 288 FBX = FAB*AX-FBB*BX 289 FBY = FAB*AY-FBB*BY 290 FBZ = FAB*AZ-FBB*BZ 291 HX(I) = HX(I)-FAX 292 HY(I) = HY(I)-FAY 293 HZ(I) = HZ(I)-FAZ 294 HX(J) = HX(J)+FAX+FBX 295 HY(J) = HY(J)+FAY+FBY 296 HZ(J) = HZ(J)+FAZ+FBZ 297 HX(K) = HX(K)-FBX 298 HY(K) = HY(K)-FBY 299 HZ(K) = HZ(K)-FBZ 300C Contributions to virials (only to projections!) 301 VIRFX = VIRFX-AX*FAX-BX*FBX 302 VIRFY = VIRFY-AY*FAY-BY*FBY 303 VIRFZ = VIRFZ-AZ*FAZ-BZ*FBZ 304 END IF 305 END DO! OF N 306CD write(*,*)TASKID,PINT*0.001*ENERF/FNOP,PINT*0.001*ENERF/FNOP 307 timea = timea+cputime(timea0) 308 RETURN 309 END 310* 311* 6. Dihedral angle contribution 312* ------------------------------- 313C 314C This subroutine calculate contribution to forces from 315C dihedral angles. Each node has own list of angles. Distribution 316C of angles over the nodes takes place in subroutine UNITS (setup.f) 317C 318C This subroutine was initially written by Andrei Komolkin 319C 320C 321C=============== GETTRS ============================================== 322C 323 SUBROUTINE GETTRS 324* 325* 6.1 front matter 326* ---------------- 327 include "prcm.h" 328 IF(NRTT.LE.0) RETURN 329 timet0 = cputime(0.d0) 330 DO M = 1,NRTT 331 TR(M) = 0.D0 332 TE(M) = 0.D0 333 END DO 334* 335* 6.2 Organize cycle over dihedral angles 336* --------------------------------------- 337 do N=1,NRRT 338 I = ITI(N) ! first atom (I) 339 ITYP = ITYPE(I) ! type 340 J = JTJ(N) ! second atom (J) 341 K = KTK(N) ! third atom (K) 342 L = LTL(N) ! forth atom (L) 343 M = ITUK(N) ! global dihedral angle number 344 NSP = NSPEC (ITYP) 345 FNSPI = 1.D0/DFLOAT(NSP) 346* 347* 6.3 Calculate dihedral angle 348* ---------------------------- 349* 6.3.1 Calculate vectors IJ,JK,KL 350 ax = sx(j)-sx(i) 351 ay = sy(j)-sy(i) 352 az = sz(j)-sz(i) 353 bx = sx(k)-sx(j) 354 by = sy(k)-sy(j) 355 bz = sz(k)-sz(j) 356 cx = sx(l)-sx(k) 357 cy = sy(l)-sy(k) 358 cz = sz(l)-sz(k) 359 call PBC(AX,AY,AZ) 360 call PBC(BX,BY,BZ) 361 call PBC(CX,CY,CZ) 362* 6.3.2 Scalar products 363 ab = ax*bx + ay*by + az*bz 364 bc = bx*cx + by*cy + bz*cz 365 ac = ax*cx + ay*cy + az*cz 366 at = ax*ax + ay*ay + az*az 367 bt = bx*bx + by*by + bz*bz 368 ct = cx*cx + cy*cy + cz*cz 369* 6.3.3 Vector products 370 axb = (at*bt)-(ab*ab) 371 bxc = (bt*ct)-(bc*bc) 372 fnum = (ab*bc)-(ac*bt) 373 den = axb*bxc 374* 6.3.4 Check that any three atoms is not on one line 375C (Otherwise contribution is zero) 376 if(den.gt.1.0d-10) then 377 den = dsqrt(den) 378C Cosine of angle: 379 co = fnum/den 380 CO = DMIN1(CO, 1.D0) 381 CO = DMAX1(CO,-1.D0) 382C Sign of angle: 383 signum = ax*(by*cz-cy*bz)+ay*(bz*cx-cz*bx)+az*(bx*cy-cx*by) 384* 6.3.5 Define angle: 385 arg = dsign(dacos(co),signum) 386 SIN1F = dsin(arg) 387 si = SIN1F 388 if( dabs(si).lt.1.0d-10 ) si = dsign(1.0d-10,si) 389* 390* 6.4 Choice from several dihedral angle types: 391* --------------------------------------------- 392* 6.4.1 Normal AMBER-type torsion angle 393 if(ITORS(M).eq.0)then 394 FCN = FT(M) 395 TEQ = RT(M) 396 MUL = NMUL(M) 397 EARG = MUL*ARG-TEQ 398 POTT = FCN*(1.D0+DCOS(EARG)) 399 DERI = -FCN*MUL*DSIN(EARG) 400* 6.4.2 MM3 force field torsional angle 401 else if(ITORS(M).eq.1)then 402 FN1 = FT1(M) 403 FN2 = FT2(M) 404 FN3 = FT3(M) 405 COS1F = dcos(ARG) 406 COS2F = 2.D0*COS1F**2-1.D0 407 COS3F = COS1F*(2.D0*COS2F-1.D0) 408 SIN2F = 2.d0*SIN1F*COS1F 409 SIN3F = SIN1F*COS2F+SIN2F*COS1F 410 POTT = FN1 * (1.D0 + COS1F) 411 X +FN2 * (1.D0 - COS2F) 412 X +FN3 * (1.D0 + COS3F) 413 DERI = -FN1*SIN1F+2.d0*FN2*SIN2F-3.d0*FN3*SIN3F 414* 6.4.3 torsional angle 415 else if(ITORS(M).eq.5)then 416 FN0 = FT(M) 417 FN1 = FT1(M) 418 FN2 = FT2(M) 419 FN3 = FT3(M) 420 FN4 = FT4(M) 421 FN5 = FT5(M) 422 TEQ = RT(M) 423 EARG = ARG-TEQ 424 COSI = DCOS(EARG) 425 COSI2 = COSI*COSI 426 COSI3 = COSI*COSI2 427 COSI4 = COSI*COSI3 428 COSI5 = COSI*COSI4 429 SINI = dsin(EARG) 430 POTT = FN0+ FN1*COSI+FN2*COSI2+ FN3*COSI3+ 431 + FN4*COSI4+ FN5*COSI5 432 DERI = -SINI*(FN1+2.d0*FN2*COSI +3.d0*FN3*COSI2+ 433 + 4.d0*FN4*COSI3+5.d0*FN5*COSI4) 434* 6.4.4 Improper dihedral angle 435 else if(ITORS(M).eq.-1)then 436 FCN = FT(M) 437 TEQ = RT(M) 438 EARG = ARG-TEQ 439 POTT = FCN*EARG**2 440 DERI = 2.d0*FCN*EARG 441 else 442 write(*,*)' Torsion type ',ITORS(M),' not supported' 443 write(*,*)' Torsion ',I,J,K,L,' on ',NAME(ITYP) 444 call FINAL 445 end if 446 TR(M) = TR(M)+abs(ARG*FNSPI) 447 TE(M) = TE(M)+POTT*FNSPI 448 PINT = PINT+POTT 449 if(IPRINT.ge.9)write(*,*)I,J,K,L,ARG*TODGR,POTT/EFACT 450* 451* 6.5 Calculate Forces: 452* --------------------- 453 if(LMOVE(ITYP))then 454 de1 = DERI/den/si 455 axb = axb/den*co 456 bxc = bxc/den*co 457* 6.5.1 X-components 458 dnum = cx*bt - bx*bc 459 dden = (ab*bx - ax*bt)*bxc 460 FFI1 = (dnum - dden) * de1 461 dnum = ((bx-ax)*bc - ab*cx ) + (2.0*ac*bx - cx*bt) 462 dden = axb*(bc*cx-bx*ct) + (ax*bt-at*bx-ab*(bx-ax))*bxc 463 FFJ1 = (dnum - dden) * de1 464 dnum = ab*bx - ax*bt 465 dden = axb*( bt*cx - bc*bx ) 466 FFL1 = (dnum - dden) * de1 467 FFK1 = -(ffi1+ffj1+ffl1) 468C Forces: 469 HX(I) = HX(I)+ffi1 470 HX(J) = HX(J)+ffj1 471 HX(K) = HX(K)+ffk1 472 HX(L) = HX(L)+ffl1 473* 6.5.2 Y-components 474 dnum = cy*bt - by*bc 475 dden = (ab*by - ay*bt)*bxc 476 FFI2 = (dnum - dden) * de1 477 dnum = ((by-ay)*bc - ab*cy ) + (2.0*ac*by - cy*bt) 478 dden = axb*(bc*cy-by*ct) + (ay*bt-at*by-ab*(by-ay))*bxc 479 FFJ2 = (dnum - dden) * de1 480 dnum = ab*by - ay*bt 481 dden = axb*( bt*cy - bc*by ) 482 FFL2 = (dnum - dden) * de1 483 FFK2 = -(ffi2+ffj2+ffl2) 484 HY(I) = HY(I)+ffi2 485 HY(J) = HY(J)+ffj2 486 HY(K) = HY(K)+ffk2 487 HY(L) = HY(L)+ffl2 488* 6.5.3 Z-components 489 dnum = cz*bt - bz*bc 490 dden = (ab*bz - az*bt)*bxc 491 FFI3 = (dnum - dden) * de1 492 dnum = ((bz-az)*bc - ab*cz ) + (2.0*ac*bz - cz*bt) 493 dden = axb*(bc*cz-bz*ct) + (az*bt-at*bz-ab*(bz-az))*bxc 494 FFJ3 = (dnum - dden) * de1 495 dnum = ab*bz - az*bt 496 dden = axb*( bt*cz - bc*bz ) 497 FFL3 = (dnum - dden) * de1 498 FFK3 = -(ffi3+ffj3+ffl3) 499 HZ(I) = HZ(I)+ffi3 500 HZ(J) = HZ(J)+ffj3 501 HZ(K) = HZ(K)+ffk3 502 HZ(L) = HZ(L)+ffl3 503* 6.5.4 Contributions to projections of virial 504 VIRFX = VIRFX-(AX+BX)*FFI1-BX*FFJ1+CX*FFL1 505 VIRFY = VIRFY-(AY+BY)*FFI2-BY*FFJ2+CY*FFL2 506 VIRFZ = VIRFZ-(AZ+BZ)*FFI3-BZ*FFJ3+CZ*FFL3 507 end if ! if(LMOVE 508 end if ! den > 0 509* 510 END DO! OF N 511 timet = timet+cputime(timet0) 512 RETURN 513 END 514* 515* 7. Special procedure for SPC water model 516* ---------------------------------------- 517C 518C This subroutine calculate intramolecular forces in 519C flexible SPC water model 520C Depending on parameter IPOT for the water, it employ 521C harmonic intramolecular potential (IPOT=1) or 522C anharmonic Morse potential (IPOT=2) 523C In the both case potential bond strength parameter defined 524C in .mol file have no meaning 525C (this is mostly because the model includes some cross-interaction 526C terms not defined in the AMBER-like potential) 527C 528C=============== STRBND =========================================== 529C 530 SUBROUTINE STRBND(ITYP) 531* 532* 7.1 Front matter 533* 534 include "prcm.h" 535C Parameters of intramolecular SPC potential (PRB 31, 2643 (1985)) 536C "parameters C, D") 537 PARAMETER (EEP = -884.63D0 , GGP = 467.31D0 ) 538C BBSA : if OH bond for water exceed BBSA, warning message is issued 539C BBMA - swithch to harmonic potential 540 PARAMETER (BBMA=1.3,BBSA=1.4 ) 541C Check that it is really water molecule 542 if(NAME(ITYP)(1:3).ne.'H2O'.and.NAME(ITYP)(1:3).ne.'h2o')then 543 write(*,*) 544 +'Potential of type 1 or 2 can be used only for water: ', 545 +' mol. name H2O*' 546 call FINAL 547 end if 548* 549* 7.2 Prepare parameters for the calculations 550* ------------------------------------------ 551C bond number for waters 552 I1 = IADB(ITYP) 553 I2 = I1+1 554 I3 = I2+1 555C Initial and last water molecules 556 IBEG = IADDR(ITYP)+1 557 IEND = IADDR(ITYP+1) 558 NSP = NSPEC(ITYP) 559 NSS = NSITS(ITYP) ! should be allways 3 560 FNSPI = 1.D0/DFLOAT(NSP) 561C Equilibrium distances 562 ROH1 = RB(I1) 563 ROH2 = RB(I2) 564 RH1H2 = RB(I3) 565C Convert potential parameters to internal units 566 FACTOR = 1.D3/AVSNO/UNITE 567* AAA = AAP*FACTOR 568* BBB = BBP*FACTOR 569* DDD = DDP 570* write(*,*)DDD,RO(I1),RO(I2) 571 IPT = IPOT(ITYP) 572 IF(IPT.EQ.2) THEN 573 AAA = DB(I1) 574 BBB = DB(I2) 575 else 576 AAA = FB(I1) 577 BBB = FB(I2) 578 END IF 579* CCC = CCP*FACTOR 580 CCC = FB(I3) 581 EEE = EEP*FACTOR 582 FFF = EEP*FACTOR 583 GGG = GGP*FACTOR 584 BBM = BBMA 585 BBS = BBSA 586 AADD = AAA*RO(I1) 587 BBDD = BBB*RO(I2) 588C These are atom numbers for the first water molecule 589 ISHF = ISADDR(ITYP) 590 IBB = ISHF+1 591 JBB = ISHF+2 592 KBB = ISHF+3 593* 594* 7.3 Case of anharmonic Morse potential 595* --------------------------------------- 596 IF(IPT.EQ.2) THEN 597* 7.3.1 Cycle over molecules 598C (each node takes care of its own set of molecules) 599 DO N = TASKID+1,NSP,NUMTASK 600C Atom numbers 601 I = (N-1)*NSS+IBB 602 J = (N-1)*NSS+JBB 603 K = (N-1)*NSS+KBB 604* 7.3.2 Bond lengths calculation 605 AX = SX(J)-SX(I) 606 AY = SY(J)-SY(I) 607 AZ = SZ(J)-SZ(I) 608 call PBC(AX,AY,AZ) 609C 610 BX = SX(K)-SX(I) 611 BY = SY(K)-SY(I) 612 BZ = SZ(K)-SZ(I) 613 call PBC(BX,BY,BZ) 614C 615 CX = SX(K)-SX(J) 616 CY = SY(K)-SY(J) 617 CZ = SZ(K)-SZ(J) 618 call PBC(CX,CY,CZ) 619C 620 RRA = AX*AX+AY*AY+AZ*AZ 621 RRB = BX*BX+BY*BY+BZ*BZ 622 RRC = CX*CX+CY*CY+CZ*CZ 623 ASS = DSQRT(RRA) 624 BSS = DSQRT(RRB) 625 CSS = DSQRT(RRC) 626 if(ASS.gt.BBS)then 627 write(*,'(a,f14.6,2i7,i10,i5)') 628 + ' !!! OH-bond length ',ASS,I,J,MSTEP,TASKID 629 write(*,*)SX(I),SY(I),SZ(I) 630 write(*,*)SX(J),SY(J),SZ(J) 631CD call XMOLDUMP 632CD stop 633 end if 634 if(BSS.gt.BBS)then 635 write(*,'(a,f14.6,2i7,i10,i5)') 636 + ' !!! OH-bond length ',BSS,I,K,MSTEP,TASKID 637 write(*,*)SX(I),SY(I),SZ(I) 638 write(*,*)SX(K),SY(K),SZ(K) 639CD call XMOLDUMP 640CD stop 641 end if 642* 7.3.3 Energy calculations 643 SA = ASS-ROH1 644 SB = BSS-ROH2 645 SH = CSS-RH1H2 646C 647 EXPA = DEXP(-RO(I1)*SA) 648 EXPAM = 1.D0-EXPA 649 EXPB = DEXP(-RO(I2)*SB) 650 EXPBM = 1.D0-EXPB 651C 652 WB1 = AAA*EXPAM**2+EEE*SA*SH 653 WB2 = BBB*EXPBM**2+FFF*SB*SH 654 WB3 = CCC*SH*SH +GGG*SA*SB 655 WB = WB1+WB2+WB3 656* 7.3.4 Fill arrays for average energy and bond length 657 BR(I1) = BR(I1)+ASS*FNSPI 658 BR(I2) = BR(I2)+BSS*FNSPI 659 BR(I3) = BR(I3)+CSS*FNSPI 660C 661 BE(I1) = BE(I1)+WB1*FNSPI 662 BE(I2) = BE(I2)+WB2*FNSPI 663 BE(I3) = BE(I3)+WB3*FNSPI 664C total intramolec. energy 665 PINT = PINT+WB1+WB2+WB3 666* 7.3.5 Force calculations 667 if(LMOVE(ITYP))then 668 if(ASS.gt.BBM)then 669C We do not have intention to dissociate water... 670 FB3A = -2.*SA*AADD*RO(I1) 671 else 672 FB3A = -2.*AADD*EXPA*EXPAM 673 end if 674 FB3A = FB3A - EEE*SH-GGG*SB 675 if(BSS.gt.BBM)then 676 FB3B = -2.*SB*BBDD*RO(I1) 677 else 678 FB3B = -2.*BBDD*EXPB*EXPBM 679 end if 680 FB3B = FB3B - FFF*SH-GGG*SA 681 FB3C = -2.0*CCC*SH - EEE*SA-FFF*SB 682* 683 FB3A = FB3A/ASS 684 FB3B = FB3B/BSS 685 FB3C = FB3C/CSS 686* 687 AXX = FB3A*AX 688 AYY = FB3A*AY 689 AZZ = FB3A*AZ 690* 691 BXX = FB3B*BX 692 BYY = FB3B*BY 693 BZZ = FB3B*BZ 694* 695 CXX = FB3C*CX 696 CYY = FB3C*CY 697 CZZ = FB3C*CZ 698* 699 HX(I) = HX(I)-AXX-BXX 700 HY(I) = HY(I)-AYY-BYY 701 HZ(I) = HZ(I)-AZZ-BZZ 702* 703 HX(J) = HX(J)+AXX -CXX 704 HY(J) = HY(J)+AYY -CYY 705 HZ(J) = HZ(J)+AZZ -CZZ 706* 707 HX(K) = HX(K) +BXX+CXX 708 HY(K) = HY(K) +BYY+CYY 709 HZ(K) = HZ(K) +BZZ+CZZ 710* 7.3.6 Virial calculations 711 VIX = AX*AXX + BX*BXX + CX*CXX 712 VIY = AY*AYY + BY*BYY + CY*CYY 713 VIZ = AZ*AZZ + BZ*BZZ + CZ*CZZ 714 VIRB = VIRB + VIX+VIY+VIZ 715 VIRFX = VIRFX+VIX 716 VIRFY = VIRFY+VIY 717 VIRFZ = VIRFZ+VIZ 718 end if 719 END DO! OF N 720* 721 ELSE 722* 723* 7.4 Harmonic flexible SPC model 724* --------------------------------- 725* 7.4.1 Cycle over water molecules 726 DO N = TASKID+1,NSP,NUMTASK 727 I = (N-1)*NSS+IBB 728 J = (N-1)*NSS+JBB 729 K = (N-1)*NSS+KBB 730* 7.4.2 Bond length calculations 731 AX = SX(J)-SX(I) 732 AY = SY(J)-SY(I) 733 AZ = SZ(J)-SZ(I) 734 BX = SX(K)-SX(I) 735 BY = SY(K)-SY(I) 736 BZ = SZ(K)-SZ(I) 737 CX = SX(K)-SX(J) 738 CY = SY(K)-SY(J) 739 CZ = SZ(K)-SZ(J) 740 RRA = AX*AX+AY*AY+AZ*AZ 741 RRB = BX*BX+BY*BY+BZ*BZ 742 RRC = CX*CX+CY*CY+CZ*CZ 743 ASS = DSQRT(RRA) 744 BSS = DSQRT(RRB) 745 CSS = DSQRT(RRC) 746* 7.4.3 Deviations 747 SA = ASS-ROH1 748 SB = BSS-ROH2 749 SH = CSS-RH1H2 750* 7.4.4 Energy calculations 751 WB1 = AAA*SA*SA+EEE*SA*SH 752 WB2 = BBB*SB*SB+FFF*SB*SH 753 WB3 = CCC*SH*SH+GGG*SA*SB 754 WB = WB1+WB2+WB3 755 BR(I1) = BR(I1)+ASS*FNSPI 756 BR(I2) = BR(I2)+BSS*FNSPI 757 BR(I3) = BR(I3)+CSS*FNSPI 758 BE(I1) = BE(I1)+WB1*FNSPI 759 BE(I2) = BE(I2)+WB2*FNSPI 760 BE(I3) = BE(I3)+WB3*FNSPI 761 PINT = PINT+WB1+WB2+WB3 762* 7.4.5 Force calculations 763 if(LMOVE(ITYP))then 764 FB3A = -2.0*AAA*SA - EEE*SH-GGG*SB 765 FB3B = -2.0*BBB*SB - FFF*SH-GGG*SA 766 FB3C = -2.0*CCC*SH - EEE*SA-FFF*SB 767 FB3A = FB3A/ASS 768 FB3B = FB3B/BSS 769 FB3C = FB3C/CSS 770 AXX = FB3A*AX 771 AYY = FB3A*AY 772 AZZ = FB3A*AZ 773 BXX = FB3B*BX 774 BYY = FB3B*BY 775 BZZ = FB3B*BZ 776 CXX = FB3C*CX 777 CYY = FB3C*CY 778 CZZ = FB3C*CZ 779 HX(I) = HX(I)-AXX-BXX 780 HY(I) = HY(I)-AYY-BYY 781 HZ(I) = HZ(I)-AZZ-BZZ 782 HX(J) = HX(J)+AXX -CXX 783 HY(J) = HY(J)+AYY -CYY 784 HZ(J) = HZ(J)+AZZ -CZZ 785 HX(K) = HX(K) +BXX+CXX 786 HY(K) = HY(K) +BYY+CYY 787 HZ(K) = HZ(K) +BZZ+CZZ 788* 7.4.6 Virial calculations 789 VIX = AX*AXX + BX*BXX + CX*CXX 790 VIY = AY*AYY + BY*BYY + CY*CYY 791 VIZ = AZ*AZZ + BZ*BZZ + CZ*CZZ 792 VIRB = VIRB + VIX+VIY+VIZ 793 VIRFX = VIRFX+VIX 794 VIRFY = VIRFY+VIY 795 VIRFZ = VIRFZ+VIZ 796 end if ! if(LMOVE 797 END DO! OF N 798 END IF 799 RETURN 800 END 801* 802* 803* 8. Calculation of LJ and electrostatic forces - small distances 804* ---------------------------------------------------------------- 805C This procedure calculates forces due to Lennard-Jones and real-space 806C electrostatics interactions for atoms on small distances: 807C closer than SHORT. These forces are calculated each short time step. 808C Interactions between atoms with distance more than SHORT are 809C calculated each short time step (see subr. FORCES, p.8) 810C 811C======================= LOCALF ========================================== 812C 813 SUBROUTINE LOCALF 814 include "prcm.h" 815 times0 = cputime(0.d0) 816 IF(MBSH.le.0) RETURN 817* 818* 8.1 Organize cycle over atom pairs 819* ----------------------------------- 820 DO INP = 1,MBSH 821 ISP0=NBS1(INP) 822 ISP=iabs(ISP0) ! atom 823 JSP=NBS2(INP) 824 ITYP=ITYPE(ISP) ! type 825 JTYP=ITYPE(JSP) 826 ISB =NSITE(ISP) ! site 827 JSB =NSITE(JSP) 828 IMOL=NNUM(ISP) ! molecule 829 JMOL=NNUM(JSP) 830 ISN =INBT(ISB) ! non-bonded type 831 JSN =INBT(JSB) 832* 833* 8.2 define interaction parameters for the given pair 834* ----------------------------------------------------- 835 MTR=MDX(ITYP,JTYP) ! type-pair index 836* same molecule 837 if(IMOL.eq.JMOL)then 838* 1-4 neighbours 839 if(ISP0.lt.0)then 840 QIJ=CHARGE(ISB)*CHARGE(JSB)*COULF*C14EL(ITYP) 841 A6=A6_14(ISN,JSN)*C14LJ(ITYP) 842 B12=B12_14(ISN,JSN)*C14LJ(ITYP) 843* other neigbours (A6LJ0 are kept constant in case of EE) 844 else 845 QIJ=CHARGE(ISB)*CHARGE(JSB)*COULF 846 A6=A6LJ0(ISN,JSN) 847 B12=B12LJ0(ISN,JSN) 848 end if 849 RAE = 0. 850* different molecules (A6LJ can be scaled for EE particles) 851 else 852 QIJ = Q(ISP)*Q(JSP)*COULF 853 A6=A6LJ(ISN,JSN) 854 B12=B12LJ(ISN,JSN) 855 if(LMDEE)then 856 if(IEE(ITYP)*IEE(JTYP).eq.0)then 857 RAE = RAEX 858 else 859 RAE = 0. 860 end if 861 end if 862 end if 863* 864* 8.3 Calculate distances between the atoms 865* ------------------------------------------ 866* 8.3.1 Vector between the atom pair 867 DX = SX(ISP)-SX(JSP) 868 DY = SY(ISP)-SY(JSP) 869 DZ = SZ(ISP)-SZ(JSP) 870* 8.3.2 Apply periodic boundary conditions 871* call PBC(DX,DY,DZ) 872 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 873 if(abs(DZ) > HBOXL) DZ = DZ - sign(BOZL,DZ) 874 if(LHEX)then 875C hexagonal periodic cell 876 XY=BOXYC*DX 877 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 878 DY=DY-BOYL 879 DX=DX-HBOXL 880 XY=BOXYC*DX 881 end if 882 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 883 DY=DY+BOYL 884 DX=DX+HBOXL 885 XY=BOXYC*DX 886 end if 887 if(DY.gt.BOXY3+XY)then 888 DY=DY-BOYL 889 DX=DX+HBOXL 890 XY=BOXYC*DX 891 end if 892 if(DY.lt.-BOXY3+XY)then 893 DY=DY+BOYL 894 DX=DX-HBOXL 895 end if 896 else 897C rectangular cell 898 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 899 end if 900 if(LOCT)then 901C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 902 CORSS = HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 903 DX=DX-sign(CORSS,DX) 904 DY=DY-sign(CORSS,DY) 905 DZ=DZ-sign(CORSS,DZ) 906 end if 907* 8.3.3 Calculate the distance and its exponents 908 RR = DX**2+DY**2+DZ**2 909 R1 = DSQRT(RR) 910 R1I = 1./R1 911 R1E = R1 + RAE 912 R1EI = 1./R1E 913* 914* 8.4 Electrostatic interaction 915* ------------------------------ 916* 8.4.1 1-4 connected pairs, or EE molecules - direct Coulomb 917 if(ISP0.le.0.or.(IMOL.eq.JMOL.and.IEE(ITYP).eq.0))then 918 EE1 = QIJ*R1I 919 if(LMOVE(ITYP))then 920 FORCE = EE1*R1I**2 921 else 922 FORCE=0. 923 end if 924 EES = EE1 925 else 926* 8.4.2 Reaction field method 927 if(LMNBL)then 928 EE1 = QIJ*(1/R1E-0.5*RFF2*R1E**2+RFF) 929 EES = EE1 930 FORCE = QIJ*(1./R1E**3+RFF2) 931C correction to virial from the reaction field is added to LJ virial 932 VIR2 = VIR2 + QIJ*(1.5*RFF2/R1E**2-RFF) 933* 8.4.3 Ewald method 934 else 935 ALPHAR = ALPHAD*R1 936 TTT = 1.D0/(1.D0+B1*ALPHAR) 937 EXP2A = DEXP(-ALPHAR**2) 938* ERFCT = erfc(TTT) 939 ERFCT = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A 940 EE1 = QIJ*R1EI 941 EES = EE1*ERFCT 942 FORCE = EE1*(TOTPI*ALPHAD*EXP2A+ERFCT*R1EI)*R1I 943 end if 944 end if 945* 8.4.4 Contributions to the electrostatic energy 946 if(IMOL.eq.JMOL)then 947 PES141(ITYP)=PES141(ITYP)+EE1 948 else 949 POTE1(MTR) = POTE1(MTR)+EE1 950 end if 951 PELS2=PELS2+EES 952* 953* 8.5 LJ interaction 954* ------------------- 955 if(B12.ne.0.d0)then 956 R2I = R1EI**2 957 R6I = R2I**3 958 ESR = ( A6+ B12*R6I)*R6I 959 FSR = (6.D0*A6+12.D0*B12*R6I)*R6I*R1EI ! |force| 960 FORCE = FORCE+FSR*R1I 961 if(IMOL.eq.JMOL)then 962 if(LMOVE(ITYP))then 963 PSR141(ITYP) = PSR141(ITYP)+ESR 964 VIR2 = VIR2+FSR*R1 965 else 966 ESR=0. 967 FORCE=0. 968 end if 969 else 970 POTL1(MTR) = POTL1(MTR)+ESR 971 VIR2 = VIR2+FSR*R1 972 end if 973 PE2 = PE2+ESR 974 end if 975* 976* 8.6 Calculate forces and virials 977* -------------------------------- 978 HX(ISP) = HX(ISP)+FORCE*DX 979 HY(ISP) = HY(ISP)+FORCE*DY 980 HZ(ISP) = HZ(ISP)+FORCE*DZ 981 HX(JSP) = HX(JSP)-FORCE*DX 982 HY(JSP) = HY(JSP)-FORCE*DY 983 HZ(JSP) = HZ(JSP)-FORCE*DZ 984 VIRFX = VIRFX+FORCE*DX**2 985 VIRFY = VIRFY+FORCE*DY**2 986 VIRFZ = VIRFZ+FORCE*DZ**2 987 END DO! OF INP 988* 989* 8.7 Intramolecular correction for "molecular" virial 990* ----------------------------------------------------- 991 DO N = 1,NSTOT 992 M = NNUM(N) 993 WIRSS = WIRSS-HX(N)*(SX(N)-X(M))-HY(N)*(SY(N)-Y(M))- 994 + HZ(N)*(SZ(N)-Z(M)) 995 END DO! OF N 996 times = times+cputime(times0) 997 RETURN 998 END 999* 1000* 9. Calculation of LJ and electrostatic forces - medium distances 1001* ---------------------------------------------------------------- 1002C This procedure calculates forces due to Lennard-Jones and real-space 1003C electrostatics interactions for atoms on "medium" distances: 1004C between SHORT and RCUT. These forces are calculated each long time step. 1005C Interactions between atoms with distance less than SHORT are 1006C calculated each short time step (see subr. LOCALF, p.8) 1007C Interactions out RCUT are assumed insignificant and do not calculated 1008C at all (long-range electrostatic forces are calculated by Ewald 1009C method, see subr. FURIER) 1010C Program flow is mostly identical to the previous section 1011C 1012C=============== SLOW FORCES ============================================= 1013C 1014 SUBROUTINE FORCES 1015 include "prcm.h" 1016 timel0 = cputime(0.d0) 1017C if nothing to do .... 1018 IF(MBLN.le.0) go to 100 1019* 1020* 9.1 Organize cycle over atom pairs 1021* ---------------------------------- 1022C List of atom pairs is determined in subroutine CHCNB 1023c This list is local for each node 1024 DO INP = 1,MBLN ! pair number 1025 ISP0=NBL1(INP) ! atom number 1026 JSP=NBL2(INP) 1027 ISP=iabs(ISP0) 1028 ITYP=ITYPE(ISP) ! type number 1029 JTYP=ITYPE(JSP) 1030 ISB =NSITE(ISP) ! site number 1031 JSB =NSITE(JSP) 1032 IMOL=NNUM(ISP) ! molecule number 1033 JMOL=NNUM(JSP) 1034 ISN =INBT(ISB) ! non-bonded type 1035 JSN =INBT(JSB) 1036* 1037* 9.2 Calculate interaction parameters for the given atom pair 1038* ------------------------------------------------------------ 1039 MTR=MDX(ITYP,JTYP) 1040* same molecule 1041 if(IMOL.eq.JMOL)then 1042* 1-4 neighbours 1043 if(ISP0.lt.0)then 1044 QIJ=CHARGE(ISB)*CHARGE(JSB)*COULF*C14EL(ITYP) 1045 A6=A6_14(ISN,JSN)*C14LJ(ITYP) 1046 B12=B12_14(ISN,JSN)*C14LJ(ITYP) 1047* other neigbours 1048 else 1049 QIJ=CHARGE(ISB)*CHARGE(JSB)*COULF 1050 A6=A6LJ0(ISN,JSN) 1051 B12=B12LJ0(ISN,JSN) 1052 end if 1053 RAE=0. 1054* different molecules 1055 else 1056 QIJ = Q(ISP)*Q(JSP)*COULF 1057 A6=A6LJ(ISN,JSN) 1058 B12=B12LJ(ISN,JSN) 1059 if(LMDEE)then 1060 if(IEE(ITYP)*IEE(JTYP).eq.0)then 1061 RAE = RAEX 1062 else 1063 RAE = 0. 1064 end if 1065 end if 1066 end if 1067* 1068* 9.3 Calculate the distance between the atoms 1069* -------------------------------------------- 1070* 9.3.1 Define the vector 1071 DX = SX(ISP)-SX(JSP) 1072 DY = SY(ISP)-SY(JSP) 1073 DZ = SZ(ISP)-SZ(JSP) 1074* 9.3.2 Apply periodic boundary conditions 1075* call PBC(DX,DY,DZ) 1076C Inline periodic boundary conditions 1077 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 1078 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 1079 if(LHEX)then 1080C hexagonal periodic cell) 1081 XY=BOXYC*DX 1082 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 1083 DY=DY-BOYL 1084 DX=DX-HBOXL 1085 XY=BOXYC*DX 1086 end if 1087 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 1088 DY=DY+BOYL 1089 DX=DX+HBOXL 1090 XY=BOXYC*DX 1091 end if 1092 if(DY.gt.BOXY3+XY)then 1093 DY=DY-BOYL 1094 DX=DX+HBOXL 1095 XY=BOXYC*DX 1096 end if 1097 if(DY.lt.-BOXY3+XY)then 1098 DY=DY+BOYL 1099 DX=DX-HBOXL 1100 end if 1101 else ! not.LHEX 1102C rectangular cell 1103 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 1104 end if 1105 if(LOCT)then 1106C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 1107 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 1108 DX=DX-sign(CORSS,DX) 1109 DY=DY-sign(CORSS,DY) 1110 DZ=DZ-sign(CORSS,DZ) 1111 end if 1112* 9.3.3 Calculate the distance 1113 RR = DX**2+DY**2+DZ**2 1114 R1 = DSQRT(RR) 1115 R1I = 1./R1 1116 R1E = R1 + RAE 1117 R1EI = 1./R1E 1118* 1119* 9.4 Electrostatic interaction 1120* --------------------------------------- 1121C EE1 is contribution to type-type energy (direct Coulomb term) 1122C EES is exact contribution to the total el-st energy 1123C FORCE is |force|/r 1124C 1125* 9.4.1 1-4 connected atoms or intramolecular EE molecules 1126C Direct Coulomb forces with unscaled charges 1127 if(ISP0.le.0.or.(IMOL.eq.JMOL.and.IEE(ITYP).eq.0))then 1128 EE1 = QIJ/R1 1129 if(LMOVE(ITYP))then 1130 FORCE = EE1/R1**2 1131 else 1132 FORCE=0. 1133 end if 1134 EES = EE1 1135 else 1136 if(LMNBL)then 1137* 9.4.2 Reaction field method 1138 EE1 = QIJ*(1/R1E-0.5*RFF2*R1E**2+RFF) 1139 EES = EE1 1140 FORCE = QIJ*(1/(R1E**2*R1)+RFF2) 1141C correction to virial from the reaction field is added to LJ virial 1142 VIR1 = VIR1 + QIJ*(1.5*RFF2*RR-RFF) 1143 else 1144* 9.4.3 Ewald method (real space term) 1145C Polynom expansion for erfc is used 1146 ALPHAR = ALPHAD*R1 1147 TTT = 1.D0/(1.D0+B1*ALPHAR) 1148 EXP2A = DEXP(-ALPHAR**2) 1149* ERFCT = erfc(TTT) 1150 ERFCT = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A 1151 EE1 = QIJ*R1EI 1152 EES = EE1*ERFCT 1153 FORCE = EE1*(TOTPI*ALPHAD*EXP2A+ERFCT*R1EI)*R1I 1154 end if 1155 end if 1156* 9.4.4 Case of EE molecule (internal electrostatic with normal charges) 1157 if(IMOL.eq.JMOL)then 1158 PES141(ITYP) = PES141(ITYP)+EE1 1159 else 1160 POTE1(MTR) = POTE1(MTR)+EE1 1161 end if 1162 PELS1=PELS1+EES 1163* 1164* 9.5 LJ interactions 1165* -------------------- 1166 if(B12.ne.0.d0)then 1167 R2I = R1EI**2 1168 R6I = R2I*R2I*R2I 1169 ESR = ( A6+ B12*R6I)*R6I 1170 FSR = (6.D0*A6+12.D0*B12*R6I)*R6I*R1EI 1171 FORCE = FORCE + FSR*R1I 1172 if(IMOL.eq.JMOL)then 1173 if(LMOVE(ITYP))VIR1 = VIR1+FSR*R1 1174 PSR141(ITYP) = PSR141(ITYP)+ESR 1175 else 1176 POTL1(MTR) = POTL1(MTR)+ESR 1177 VIR1 = VIR1+FSR*R1 1178 end if 1179 PE1 = PE1 + ESR 1180 end if 1181* 1182* 9.6 Contributions to virial projections 1183* ----------------------------------------- 1184 if(IMOL.ne.JMOL.or.LMOVE(ITYP))then 1185 VIRX = VIRX+FORCE*DX**2 1186 VIRY = VIRY+FORCE*DY**2 1187 VIRZ = VIRZ+FORCE*DZ**2 1188 end if 1189CD write(*,*)ISP,JSP,IMOL,JMOL,R1,ESR*PERMOL,EES*PERMOL 1190* 1191* 9.7 Calculate forces 1192* --------------------- 1193 GX(ISP) = GX(ISP)+FORCE*DX 1194 GY(ISP) = GY(ISP)+FORCE*DY 1195 GZ(ISP) = GZ(ISP)+FORCE*DZ 1196 GX(JSP) = GX(JSP)-FORCE*DX 1197 GY(JSP) = GY(JSP)-FORCE*DY 1198 GZ(JSP) = GZ(JSP)-FORCE*DZ 1199 END DO! OF INP 1200 100 timel = timel+cputime(timel0) 1201 RETURN 1202 END 1203* 1204*============================================================= 1205* 1206* 10. Ewald method. Reciprocal space contribution 1207* ----------------------------------------------- 1208C This procedure calculates contribution to the energy and forces 1209C from long-range part of electrostatic interactions, calculated 1210C in the reciprocal space by the Furie transform 1211C 1212C Subroutines: 1213C FURIR - the main procedure for Ewald sum calculation 1214C INIFUR - initialisation of arrays for Ewald sum calculation 1215C=============== FURIR ======================================= 1216C 1217 SUBROUTINE FURIR 1218* 1219* 10.1 Definitions 1220* ----------------- 1221 include "prcm.h" 1222 parameter (FKTX=1.05,FKTS=0.95) 1223 parameter(NPUSH=10000) 1224 DIMENSION TSIN(NTOT),TCOS(NTOT) 1225 data INIT/0/ 1226* 1227* 10.2 Initialization 1228* -------------------- 1229 timef0 = cputime(0.d0) 1230 if(INIT.eq.0)then 1231 INIT=1 1232 NFPUSH=NPUSH/NSTOT+1 1233* 10.2.1 remember box size 1234 BOXLO=BOXL 1235 BOYLO=BOYL 1236 BOZLO=BOZL 1237* 10.2.2 Initialize arrays for the calculations 1238 if(ALPHAD.gt.0.d0)then 1239 call INIFUR 1240 ALP2=ALPHAD**2 1241 if(TASKID.eq.MAST.and.IPRINT.ge.4)write(*,*) 1242 +'Ewald sum: Num. of k-vectors ',NKV 1243* 10.2.3 Case when Ewald method is not used 1244 else 1245 NKV=0 1246 if(TASKID.eq.MAST.and..not.LMNBL)write(*,*) 1247 +'No Ewald' 1248 end if 1249 end if 1250* 10.2.4 If box size changed too strong, reinitialize 1251 DBX=BOXL/BOXLO 1252 DBY=BOYL/BOYLO 1253 DBZ=BOZL/BOZLO 1254 if(DBX.gt.FKTX.or.DBX.lt.FKTS.or.DBY.gt.FKTX.or. 1255 * DBY.lt.FKTS.or.DBZ.gt.FKTX.or.DBZ.lt.FKTS)then 1256 BOXLO=BOXL 1257 BOYLO=BOYL 1258 BOZLO=BOZL 1259 call INIFUR 1260 ALP2=ALPHAD**2 1261 if(TASKID.eq.MAST.and.IPRINT.ge.4)write(*,*) 1262 +'Box size changed. New num. of k-vectors ',NKV 1263C for expanded ensemble 1264 if(LMDEE)call DIFF_LRCOR 1265 end if 1266 IF(NKV.eq.0.or.NTYPES.EQ.1.AND.NSPEC(1).EQ.1) RETURN 1267* 1268* 10.3 Reciprocal space Ewald sum 1269* -------------------------------- 1270 QPE = 0.D0 1271 CFUR=4.d0*PI*COULF/VOL 1272* 10.3.1 Start cycle over vectors in the reciprocal space 1273 IBEG=NUMTASK-TASKID 1274 if(LVISUAL)then 1275 if(mod(IKV,NFPUSH).eq.0)write(*,'(a)')'@mm ' 1276 end if 1277 do IKV=IBEG,NKV,NUMTASK 1278C These are coordinates of reciprocal vectors 1279 RXV=RKX(IKV)/BOXL 1280 RYV=RKY(IKV)/BOYL 1281 RZV=RKZ(IKV)/BOZL 1282 RK2=RXV**2+RYV**2+RZV**2 1283 SCS=0. 1284 SSN=0. 1285* 10.3.2 Calculate Sin(rk) and Cos(rk) and their sums for all the atoms 1286 do I=1,NSTOT 1287C Scalar product r.k 1288 SCP=RXV*SX(I)+RYV*SY(I)+RZV*SZ(I) 1289C sin and cos are taken from look-up tables 1290 QQ=Q(I) 1291 SINSC=QQ*sin(SCP) 1292 COSSC=QQ*cos(SCP) 1293C sin and cos calculated from look-up tables 1294C if(SCP.le.0.d0)then 1295C XT=-SCP*DTAB 1296C IID=1 1297C else 1298C XT=SCP*DTAB 1299C IID=0 1300C end if 1301C IOT=idint(XT) 1302C IU=mod(IOT,LTAB) 1303C IU1=IU+1 1304C DD=XT-dfloat(IOT) 1305C DD1=1.d0-DD 1306C SINSC=Q(I)*(DD1*SINT(IU)+DD*SINT(IU1)) 1307C COSSC=Q(I)*(DD1*COST(IU)+DD*COST(IU1)) 1308C if(IID.ne.0)SINSC=-SINSC 1309* 10.3.3 remember sin and cos in temporarly arrays (for forces) 1310 TSIN(I)=SINSC 1311 TCOS(I)=COSSC 1312C These are sums 1313 SSN=SSN+SINSC 1314 SCS=SCS+COSSC 1315 end do 1316* 10.3.4 Calculate energy 1317 AKC=CFUR*exp(-0.25d0*RK2/ALP2)/RK2 1318 QPE=QPE+AKC*(SSN**2+SCS**2) 1319* 10.3.5 Calculate contribution to forces 1320 AK2=2.d0*AKC 1321 do I=1,NSTOT 1322 QFOR=AK2*(SCS*TSIN(I)-SSN*TCOS(I)) 1323 GX(I)=GX(I)+RXV*QFOR 1324 GY(I)=GY(I)+RYV*QFOR 1325 GZ(I)=GZ(I)+RZV*QFOR 1326 end do 1327* 10.3.6 Calculate contribution to virial projections 1328 AKV=0.5*AKC*(4.*ALP2+RK2)/(ALP2*RK2) 1329 VIRX=VIRX-RXV**2*AKV*(SSN**2+SCS**2) 1330 VIRY=VIRY-RYV**2*AKV*(SSN**2+SCS**2) 1331 VIRZ=VIRZ-RZV**2*AKV*(SSN**2+SCS**2) 1332 end do 1333 VIRX=VIRX+QPE 1334 VIRY=VIRY+QPE 1335 VIRZ=VIRZ+QPE 1336 PELS1=PELS1+QPE 1337 timef = timef+cputime(timef0) 1338 return 1339 end 1340C================================================================= 1341* 1342* 11. Initialisation of Ewald sum calculations 1343* -------------------------------------------- 1344* 1345C This subroutine defines a set of vectors in the reciprocal space 1346C which will be used for Ewald sum calculations in dependence of geometry 1347C Basis vectors in the reciprocal space are determined from vector 1348C products: Kx = Ly x Lz and so on (Lxyz are basis vectros in real space) 1349C For hexagonal cell: Kx=(1,1/sqrt(3),0), Ky=(0,2/sqrt(3),0), Kz=(0,0,1) 1350C For octahedron cell: Kx=(1,0,-1), Ky=(0,1,-1), Kz=(0,0,2) 1351C (multiplied by 2*Pi/Lbox in each direction; note that for hexagonal cell 1352C BOYL defined as BOYL=sqrt(3)*Ly/2) 1353* 1354*================= INIFUR ========================================= 1355* 1356 subroutine INIFUR 1357 include "prcm.h" 1358* 1359* 11.1 Calculate look-up tables for sin and cos 1360* ---------------------------------------------- 1361C call SINTAB 1362* 1363* 11.2 Calculate parameters which set up cut-off in the reciprocal space 1364* ----------------------------------------------------------------------- 1365 ALP2=ALPHAD**2 1366 CUTK2=4.d0*ALP2*FEXP 1367 CUTK=sqrt(CUTK2) 1368 TWOPI=2.d0*PI 1369 IKV=0 1370 KMAX=CUTK*BOXL/TWOPI+1 1371 KMAY=CUTK*BOYL/TWOPI+1 1372 KMAZ=CUTK*BOZL/TWOPI+1 1373* 1374* 11.3 Calculate set of reciprocal vectors 1375* ----------------------------------------- 1376* 11.3.1 Cycle over reciprocal vectors in a cube 1377 do IKZ=-KMAZ,KMAZ 1378 do IKX=0,KMAX 1379 if(IKX.eq.0)then 1380 if(IKZ.gt.0)then 1381 KY1=0 1382 else 1383 KY1=1 1384 end if 1385 KY2=KMAY 1386 else 1387 KY1=-KMAY 1388 KY2=KMAY 1389 end if 1390 do IKY=KY1,KY2 1391* 11.3.2 Calculate cartesian coordinates of reciprocal vectors 1392 if(ICELL.eq.1)then 1393C Truncated octahedron 1394C These are cartesian coordinates of reciprocal vectors multiplied by 1395C the box size (which may change in NPT-dynamics) 1396 RXV = TWOPI*IKX 1397 RYV = TWOPI*IKY 1398 RZV = TWOPI*(-IKX+IKY+2*IKZ) 1399 else if(ICELL.eq.2)then 1400C Hexagonal 1401 RXV = TWOPI*IKX 1402 RYV = TWOPI*IKY+PI*IKX 1403 RZV = TWOPI*IKZ 1404 else ! ICELL=0 1405C Rectangular cell 1406 RXV = TWOPI*IKX 1407 RYV = TWOPI*IKY 1408 RZV = TWOPI*IKZ 1409 end if 1410* 11.3.3 Remember reciprocal vectors which are within cutoff 1411 RK2 = (RXV/BOXL)**2+(RYV/BOYL)**2+(RZV/BOZL)**2 1412 if(RK2.le.CUTK2)then 1413 IKV=IKV+1 1414 if(IKV.gt.NKVM)then 1415 write(*,*)'max number of recirpocal vectors exceeded' 1416 write(*,*)'check Ewald summation parameters or ' 1417 write(*,*)'increase NKVM in forces.f/FURIER and INIFUR' 1418 stop 'in INIFUR' 1419 end if 1420 RKX(IKV)=RXV 1421 RKY(IKV)=RYV 1422 RKZ(IKV)=RZV 1423 end if 1424 end do 1425 end do 1426 end do 1427 NKV=IKV 1428* 11.4 Ewald self-energy 1429 ABC = ALPHAD/DSQRT(PI) 1430 SPEI=0. 1431 do ITYP=1,NTYPES 1432 SPET=0. 1433C EE scaling factor for ITYP 1434 if(IEE(ITYP).eq.0)then 1435 ESCI=EC(ME) 1436 else 1437 ESCI=1.d0 1438 end if 1439 do I=1,NSITS(ITYP) 1440 IS=ISADR(ITYP)+I 1441 QII=COULF*(ESCI*CHARGE(IS))**2 1442 SPET = SPET-NSPEC(ITYP)*ABC*QII 1443 end do 1444 SPET=SPET/NUMTASK 1445 SPEEI(ITYP)=SPET 1446 SPEI=SPEI+SPET 1447 end do 1448 return 1449 end 1450* 1451* 11.4 Create look-up tables for sin and cos 1452* ------------------------------------------- 1453*=========================== SINTAB ================================= 1454* 1455C subroutine SINTAB 1456C implicit real*8 (A-H,O-Z) 1457C parameter (PI= 3.1415926535897932d0, PI2=2.d0*PI, LTAB=50000) 1458C common /TSNCS/ DTAB,SINT(0:LTAB),COST(0:LTAB) 1459C do I=0,LTAB 1460C X=dfloat(I)*PI2/LTAB 1461C sint(I)=sin(X) 1462C cost(I)=cos(X) 1463C end do 1464C DTAB=dfloat(LTAB)/PI2 1465C return 1466C end 1467C======================================================================== 1468* 1469* 12. Calculation of self-interaction term in the Ewald sum 1470* ---------------------------------------------------------- 1471C 1472C Really this procedure calculates two terms: 1473C 1) Standard self-interaction Ewald term, giving only constant 1474C contribution to the energy 1475C 2) "Molecular correction" which is due to the fact, that bound 1476C atoms in a molecule do not interact electrostatically, but 1477C reciprocal Ewald term includes these interactions also. This correction 1478C exclude contribution of reciprocal space Ewald between certain atom pairs 1479* 1480*============= ETERMS ================================================= 1481* 1482 SUBROUTINE ETERMS 1483* 12.1 Front matter 1484 include"prcm.h" 1485 data INIT/0/ 1486 timee0 = cputime(0.d0) 1487 if(LMNBL)return 1488* 1489* 12.2 Electrostatic self-interaction 1490* ------------------------------------ 1491* (constant terms SPEEI and SPEI are computed separately) 1492 do I=1,NTYPES 1493 SPEE(I)=SPEEI(I) 1494 end do 1495 SPE=SPEI 1496* 1497* 12.3 Calculate intramolecular correction 1498* ----------------------------------------- 1499* 12.3.1 Start cycle over bound atoms 1500C list of bound atoms is prepeared in BNDLST (setup.f) 1501 IBEG=NUMTASK-TASKID 1502 do ISP=IBEG,NSTOT,NUMTASK 1503 ITYP = ITYPE(ISP) 1504* this is normal case (not EE molecules) 1505 if(IEE(ITYP).ne.0)then 1506 NSP = NNBB(ISP) 1507 do J=1,NSP 1508 JSP= INBB(ISP,J) 1509 if(JSP.le.0)JSP=-JSP 1510 if(ISP.gt.JSP)then ! avoid double count! 1511 IS = NSITE(ISP) 1512 JS = NSITE(JSP) 1513 QIJ=COULF*Q(ISP)*Q(JSP) 1514 DX = SX(ISP)-SX(JSP) 1515 DY = SY(ISP)-SY(JSP) 1516 DZ = SZ(ISP)-SZ(JSP) 1517* 12.3.2 Apply periodic boundary conditions 1518 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 1519 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 1520 if(LHEX)then 1521C hexagonal periodic cell) 1522 XY=BOXYC*DX 1523 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 1524 DY=DY-BOYL 1525 DX=DX-HBOXL 1526 XY=BOXYC*DX 1527 end if 1528 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 1529 DY=DY+BOYL 1530 DX=DX+HBOXL 1531 XY=BOXYC*DX 1532 end if 1533 if(DY.gt.BOXY3+XY)then 1534 DY=DY-BOYL 1535 DX=DX+HBOXL 1536 XY=BOXYC*DX 1537 end if 1538 if(DY.lt.-BOXY3+XY)then 1539 DY=DY+BOYL 1540 DX=DX-HBOXL 1541 end if 1542 else 1543C rectangular cell 1544 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 1545 end if 1546 if(LOCT)then 1547 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 1548 DX=DX-sign(CORSS,DX) 1549 DY=DY-sign(CORSS,DY) 1550 DZ=DZ-sign(CORSS,DZ) 1551 end if 1552* call PBC(DX,DY,DZ) 1553* 12.3.3 Calculate contributions to energy and forces 1554 R2 = DX**2+DY**2+DZ**2 1555 R2I = 1.D0/R2 1556 R1 = DSQRT(R2) 1557 R1I = R1*R2I 1558 ALPHAR = ALPHAD*R1 1559* call RERFC(ALPHAR,ERFC,EXP2A) 1560 TTT = 1.D0/(1.D0+B1*ALPHAR) 1561 EXP2A = DEXP(-ALPHAR**2) 1562 ERFC = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A 1563 EES = QIJ*(ERFC-1.d0)*R1I 1564 SPE = SPE+EES 1565 SPEE(ITYP) = SPEE(ITYP)+EES 1566 if(LMOVE(ITYP))then 1567 FES = QIJ*((TOTPI*ALPHAR*EXP2A+ERFC)-1.d0)*R1I*R2I 1568 GX(ISP) = GX(ISP)+FES*DX 1569 GY(ISP) = GY(ISP)+FES*DY 1570 GZ(ISP) = GZ(ISP)+FES*DZ 1571 GX(JSP) = GX(JSP)-FES*DX 1572 GY(JSP) = GY(JSP)-FES*DY 1573 GZ(JSP) = GZ(JSP)-FES*DZ 1574 VIRX = VIRX+FES*DX**2 1575 VIRY = VIRY+FES*DY**2 1576 VIRZ = VIRZ+FES*DZ**2 1577 end if 1578 end if 1579 end do 1580 end if 1581 END DO! OF I 1582* 1583* 12.4 Calculate intramolecular correction for EE molecules 1584* ---------------------------------------------------------- 1585* 12.4.1 find EE pairs 1586* EE molecules: direct nonscaled Coulomb between all atom pairs (in parts 8 and 9) 1587* here we remove reciprocal space contribution inside EE molecules 1588* (made in FOURIER with scaled charges) 1589* find EE molecules 1590 do ITYP=1,NTYPES 1591 if(IEE(ITYP).eq.0)then 1592* cycle over molecules of this type 1593 do IM=TASKID,NSPEC(ITYP)-1,NUMTASK 1594 NSP=NSITS(ITYP) ! atoms in the molecule 1595 ISPI=ISADDR(ITYP)+NSP*IM ! first atom in mol - 1 1596 if(NSP.gt.1)then 1597 do ISM=2,NSP 1598 do JSM=1,ISM-1 1599 ISP=ISPI+ISM 1600 JSP=ISPI+JSM 1601 IS = NSITE(ISP) 1602 JS = NSITE(JSP) 1603 QIJ=COULF*Q(ISP)*Q(JSP) ! scaled charges 1604 DX = SX(ISP)-SX(JSP) 1605 DY = SY(ISP)-SY(JSP) 1606 DZ = SZ(ISP)-SZ(JSP) 1607* 12.4.2 Apply periodic boundary conditions 1608 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 1609 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 1610 if(LHEX)then 1611C hexagonal periodic cell) 1612 XY=BOXYC*DX 1613 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 1614 DY=DY-BOYL 1615 DX=DX-HBOXL 1616 XY=BOXYC*DX 1617 end if 1618 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 1619 DY=DY+BOYL 1620 DX=DX+HBOXL 1621 XY=BOXYC*DX 1622 end if 1623 if(DY.gt.BOXY3+XY)then 1624 DY=DY-BOYL 1625 DX=DX+HBOXL 1626 XY=BOXYC*DX 1627 end if 1628 if(DY.lt.-BOXY3+XY)then 1629 DY=DY+BOYL 1630 DX=DX-HBOXL 1631 end if 1632 else 1633C rectangular cell 1634 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 1635 end if 1636 if(LOCT)then 1637 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 1638 DX=DX-sign(CORSS,DX) 1639 DY=DY-sign(CORSS,DY) 1640 DZ=DZ-sign(CORSS,DZ) 1641 end if 1642* call PBC(DX,DY,DZ) 1643* 12.4.3 Calculate contributions to energy and forces 1644 R2 = DX**2+DY**2+DZ**2 1645 R2I = 1.D0/R2 1646 R1 = DSQRT(R2) 1647 R1I = R1*R2I 1648 ALPHAR = ALPHAD*R1 1649* call RERFC(ALPHAR,ERFC,EXP2A) 1650 TTT = 1.D0/(1.D0+B1*ALPHAR) 1651 EXP2A = DEXP(-ALPHAR**2) 1652 ERFC = ((((A5*TTT+A4)*TTT+A3)*TTT+A2)*TTT+A1)*TTT*EXP2A 1653 EES = QIJ*(ERFC-1.d0)*R1I 1654 SPE = SPE+EES 1655 SPEE(ITYP) = SPEE(ITYP)+EES 1656 if(LMOVE(ITYP))then 1657 FES = QIJ*((TOTPI*ALPHAR*EXP2A+ERFC)-1.d0)*R1I*R2I 1658 GX(ISP) = GX(ISP)+FES*DX 1659 GY(ISP) = GY(ISP)+FES*DY 1660 GZ(ISP) = GZ(ISP)+FES*DZ 1661 GX(JSP) = GX(JSP)-FES*DX 1662 GY(JSP) = GY(JSP)-FES*DY 1663 GZ(JSP) = GZ(JSP)-FES*DZ 1664 VIRX = VIRX+FES*DX**2 1665 VIRY = VIRY+FES*DY**2 1666 VIRZ = VIRZ+FES*DZ**2 1667 end if 1668 end do 1669 end do 1670 end if ! if(NSP... 1671 end do 1672 end if 1673 END DO! OF ITYP 1674 timee = timee+cputime(timee0) 1675 RETURN 1676 END 1677*===================================================================== 1678* 1679* 13. Calculate out-cut-off corrections to energy and pressure from 1680* LJ interactions 1681* 1682*============= ELRCLJ ================================================= 1683* 1684 SUBROUTINE ELRCLJ 1685 include"prcm.h" 1686 dimension PELRC0(NTPP),VRLRC0(NTPP) 1687 data INIT/0/ 1688 if(INIT.eq.0)then 1689 DO I = 1,MOLINT 1690 PELRC0(I) = 0.D0 1691 VRLRC0(I) = 0.D0 1692 END DO! OF I 1693 RC0 = DSQRT(RSSQ) 1694C Cycle over all site pairs 1695 DO 400 ITYP = 1 ,NTYPES 1696 ISBEG = ISADR(ITYP)+1 1697 ISEND = ISADR(ITYP +1) 1698 NSPI = NSPEC(ITYP) 1699 DO 400 IS = ISBEG,ISEND 1700 IST = INBT(IS) 1701 DO 300 JTYP = ITYP ,NTYPES 1702 if(IEE(ITYP)*IEE(JTYP).eq.0)then 1703 RC = RC0 + RAEX 1704 else 1705 RC = RC0 1706 end if 1707 MT = MDX (ITYP,JTYP) 1708 JSBEG = ISADR(JTYP)+1 1709 JSEND = ISADR(JTYP +1) 1710 NSPJ = NSPEC(JTYP) 1711 FNOPIJ = DFLOAT(NSPI*NSPJ) 1712 if(ITYP.eq.JTYP)FNOPIJ=0.5*FNOPIJ 1713 DO 300 JS = JSBEG,JSEND 1714 JST = INBT(JS) 1715 if(B12LJ(IST,JST).gt.1.d-50.and. 1716 + A6LJ(IST,JST).lt.-1.d-50)then 1717 EP = A6LJ(IST,JST)**2/B12LJ(IST,JST) 1718 SI3 = sqrt(-B12LJ(IST,JST)/A6LJ(IST,JST)) ! SI**3 1719 PELRC0(MT) = PELRC0(MT)+EP*SI3*(( 1.D0/9.D0)*(SI3/RC**3)**3 1720 X - ( 1.D0/3.D0)*(SI3/RC**3))*FNOPIJ 1721 VRLRC0(MT) = VRLRC0(MT)+EP*SI3*((-4.D0/3.D0)*(SI3/RC**3)**3 1722 X + ( 2.D0)*(SI3/RC**3))*FNOPIJ 1723 end if 1724 300 CONTINUE 1725 400 CONTINUE 1726 end if 1727 FCV = UNITP/(3.*VOL) 1728 VIRD = 0. 1729 EPD = 0. 1730 DO I = 1,MOLINT 1731 PELRC(I) = 4.*PI*PELRC0(I)/VOL 1732 VRLRC(I) = -4.*PI*VRLRC0(I)/VOL 1733 EPD = EPD+PELRC(I) 1734 VIRD=VIRD+VRLRC(I) 1735 if(INIT.eq.0.and.IPRINT.ge.5)then 1736 PRINT "('*** PE LRC ',I2,': ',F12.8,' kJ/mol')", 1737 + I,PELRC(I)*PERMOL 1738 PRINT "('*** PR LRC ',I2,': ',F12.4,' atm ')",I,VRLRC(I)*FCV 1739 INIT=1 1740 end if 1741 END DO 1742* 1743 RETURN 1744 END 1745C=================================================================== 1746* 1747* 14. Harmonic force aroun fixed position 1748* --------------------------------------- 1749C 1750C Each atom may be put in a harmonic well around 1751C some position defined in "filref" file (if LRR=.true) 1752*====================== PULLBACK ==================================== 1753* 1754 subroutine PULLBACK 1755 include "prcm.h" 1756* Drag force applied to some molecule types 1757 if(LDRAG)then 1758 do I=NAB(TASKID),NAE(TASKID) 1759 ITYP=ITYPE(I) 1760 if(LDRAGT(ITYP))then 1761 IS = NSITE(I) 1762 GZ(I)=GZ(I)+MASSD(IS)*FDRAG(ITYP) 1763 end if 1764 end do 1765 end if 1766 if((.not.LRR.or.NPULL.eq.0).and.IEXT.ne.1)return 1767 if(LRR)then 1768 ELINK=0. 1769 do IL=1,NPULL ! num of linkage point 1770 I = INPL(IL) ! num of atom 1771 if(I.le.0)write(*,*)' wrong index in PULLBACK!!!' 1772 DX=SX(I)-XPL(IL) 1773 DY=SY(I)-YPL(IL) 1774 DZ=SZ(I)-ZPL(IL) 1775C Normally, PBC are not used here 1776C call PBC(DX,DY,DZ) 1777 RR2=DX**2+DY**2+DZ**2 1778 ELINK=ELINK+RR2 1779 GX(I)=GX(I)-RLR*DX 1780 GY(I)=GY(I)-RLR*DY 1781 GZ(I)=GZ(I)-RLR*DZ 1782 end do 1783 ELNK=ELINK*PERMOL*RLR 1784 if(IPRINT.ge.6)write(*,*) 1785 +'Deviation from correct structure ',sqrt(ELINK/NPULL),ELNK 1786 end if 1787* External electrostatic field 1788 if(IEXT.eq.1)then 1789 FULTIM=TIM+NSTEP*DT 1790 EFIELD=EXAMPL*cos(EXFREQ*FULTIM*0.5/PI) 1791 do I=NAB(TASKID),NAE(TASKID) 1792 IS=NSITE(I) 1793 GZ(I)=GZ(I)+Q(I)*EFIELD 1794 end do 1795 end if 1796 return 1797 end 1798C 1799C=========================================================================== 1800* 1801* 15. Recalculation of list of neighbours 1802* --------------------------------------- 1803C This procedure calculate lists of closest and far neighbours 1804C for subroutines LOCALF and FORCES 1805C These lists must be recalculated after some time 1806C Each node calculate goes through its own list of atom pairs 1807C and creates own lists, which are used in subroutines for force calculation 1808C 1809C================ CHCNB ======================================================= 1810C 1811 subroutine CHCNB(LRDF) 1812* 15.1 Front matter 1813 include "prcm.h" 1814 parameter (RMX2=8.**2,MMOL=5) 1815 data INIT/0/ 1816 timev0 = cputime(0.d0) 1817 MBLN=0 1818 MBSH=0 1819 IOVER=0 1820* 1821* 15.2 Organize cycle over ALL atom pairs 1822* ---------------------------------------- 1823 IBEG=NUMTASK-TASKID 1824C cycle over I,J is organized is follows: 1825C it takes I>J if I and J are of the same parity and I<J if I and J have 1826C different parity. This allows nearly uniformly distribute atom pairs 1827C between the nodes 1828 do I=IBEG,NSTOT,NUMTASK 1829 IMOL=NNUM(I) 1830 ITYP=ITYPE(I) 1831* 15.2.1 Even strips 1832 do J=I-2,1,-2 1833 JMOL=NNUM(J) 1834 DX=SX(I)-SX(J) 1835 DY=SY(I)-SY(J) 1836 DZ=SZ(I)-SZ(J) 1837* 15.2.1.1 Periodic boundary conditions 1838C call PBC (DX,DY,DZ) 1839 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 1840 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 1841 if(LHEX)then 1842C hexagonal periodic cell) 1843 XY=BOXYC*DX 1844 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 1845 DY=DY-BOYL 1846 DX=DX-HBOXL 1847 XY=BOXYC*DX 1848 end if 1849 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 1850 DY=DY+BOYL 1851 DX=DX+HBOXL 1852 XY=BOXYC*DX 1853 end if 1854 if(DY.gt.BOXY3+XY)then 1855 DY=DY-BOYL 1856 DX=DX+HBOXL 1857 XY=BOXYC*DX 1858 end if 1859 if(DY.lt.-BOXY3+XY)then 1860 DY=DY+BOYL 1861 DX=DX-HBOXL 1862 end if 1863 else 1864C rectangular cell 1865 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 1866 end if 1867 if(LOCT)then 1868C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 1869 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 1870 DX=DX-sign(CORSS,DX) 1871 DY=DY-sign(CORSS,DY) 1872 DZ=DZ-sign(CORSS,DZ) 1873 end if 1874 RR2=DX**2+DY**2+DZ**2 1875* 15.2.1.2 Case of "molecular" cutoff - calculate COM distance 1876C If we do not use Ewald summation (LMNBL=.true.), 1877C it is important that at least 1878C small molecules (less than MMOL atoms) where inside or outside 1879C cut-off radius as a whole 1880 if(LMNBL)then 1881 NMI=NSITS(ITYP) 1882 NMJ=NSITS(ITYPE(J)) 1883 if(NMI.gt.MMOL.and.NMJ.gt.MMOL)then 1884 RRM=RR2 1885 go to 110 1886 else if(NMI.gt.MMOL.and.NMJ.le.MMOL)then 1887 DX=SX(I)-X(JMOL) 1888 DY=SY(I)-Y(JMOL) 1889 DZ=SZ(I)-Z(JMOL) 1890 else if(NMJ.gt.MMOL)then 1891 DX=X(IMOL)-SX(J) 1892 DY=Y(IMOL)-SY(J) 1893 DZ=Z(IMOL)-SZ(J) 1894 else 1895 DX=X(IMOL)-X(JMOL) 1896 DY=Y(IMOL)-Y(JMOL) 1897 DZ=Z(IMOL)-Z(JMOL) 1898 end if 1899* call PBC (DX,DY,DZ) 1900 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 1901 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 1902 if(LHEX)then 1903C hexagonal periodic cell) 1904 XY=BOXYC*DX 1905 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 1906 DY=DY-BOYL 1907 DX=DX-HBOXL 1908 XY=BOXYC*DX 1909 end if 1910 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 1911 DY=DY+BOYL 1912 DX=DX+HBOXL 1913 XY=BOXYC*DX 1914 end if 1915 if(DY.gt.BOXY3+XY)then 1916 DY=DY-BOYL 1917 DX=DX+HBOXL 1918 XY=BOXYC*DX 1919 end if 1920 if(DY.lt.-BOXY3+XY)then 1921 DY=DY+BOYL 1922 DX=DX-HBOXL 1923 end if 1924 else 1925C rectangular cell 1926 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 1927 end if 1928 if(LOCT)then 1929C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 1930 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 1931 DX=DX-sign(CORSS,DX) 1932 DY=DY-sign(CORSS,DY) 1933 DZ=DZ-sign(CORSS,DZ) 1934 end if 1935C square distance betwen COM of small molecules 1936 RRM=DX**2+DY**2+DZ**2 1937 else ! .not.LMNBL 1938 RRM=RR2 1939 end if ! if(LMNBL 1940* 15.2.1.3 Prepare list of intramolecular non-bonded pairs 1941 110 if(IMOL.eq.JMOL)then ! bound atoms should not be in the list 1942 if(RR2.le.RMX2)then 1943C Check that the atoms are non-bound. Check using binding list of the 1944C atom which has less bound atoms 1945 do K=1,NNBB(I) 1946C 1-3 bound: not in the list 1947 if(INBB(I,K).eq.J)go to 125 1948 if(INBB(I,K).eq.-J)then 1949C 1-4 bound: put into list with negative first index 1950C (these atom pairs may have a special way of force calculations) 1951 MBSH=MBSH+1 1952 if(MBSH.gt.NBSMAX)write(*,*) 1953 + '!!! list of closest neigbours overfilled. atoms ',I,J 1954 NBS1(MBSH)=-I 1955 NBS2(MBSH)=J 1956 go to 125 1957 end if 1958 end do 1959 end if ! (RR2.le.RMAX2 1960 end if ! (IMOL.eq.JMOL 1961* 15.2.1.4 Pick up non-bound atoms to the list of neigbours 1962 if(RRM.le.SHORT2)then 1963 MBSH=MBSH+1 1964 if(MBSH.gt.NBSMAX)write(*,*) 1965 + '!!! list of closest neigbours overfilled. atoms ',I,J 1966 NBS1(MBSH)=I 1967 NBS2(MBSH)=J 1968 else if(RRM.le.RSSQ)then 1969 MBLN=MBLN+1 1970 if(MBLN.gt.NBLMAX)then 1971 if(IOVER.eq.0)write(*,*) 1972 + '!! list of far neigbours overfilled' 1973 IOVER=1 1974 end if 1975 if(IOVER.eq.0)then 1976 NBL1(MBLN)=I 1977 NBL2(MBLN)=J 1978 end if 1979 end if 1980* 15.2.1.5 RDF calculation 1981C Only at this point we go through all the atom pairs 1982 125 continue 1983*X IF(LGRC.and.RR2.le.RDFCUT2) THEN !=====> LGR 1984*X ISB = NSITE(J) 1985*X JSB = NSITE(I) 1986*XC find RDF pair 1987*X if(NGRI(ISB).le.NGRI(JSB))then 1988*X do JS=1,NGRI(ISB) 1989*X if(IGRI(ISB,JS).eq.JSB)then 1990*X IP = MGRI(ISB,JS) 1991*X go to 127 1992*X end if 1993*X end do 1994*X else 1995*X do JS=1,NGRI(JSB) 1996*X if(IGRI(JSB,JS).eq.ISB)then 1997*X IP = MGRI(JSB,JS) 1998*X go to 127 1999*X end if 2000*X end do 2001*X end if 2002*X go to 128 2003*X 127 R1 = sqrt(RR2) 2004*X GG = R1*DRDFI 2005*X MM = IDINT(GG)+1 2006*X if(MM.le.MAX)IRDF(MM,IP) = IRDF(MM,IP)+1 2007*X end if ! of LGRC 2008 128 continue 2009 end do 2010* 15.2.2 Odd strips 2011C Everythig here is a complete analogue to p 15.2.1 2012 do J=I+1,NSTOT,2 2013 JMOL=NNUM(J) 2014 DX=SX(I)-SX(J) 2015 DY=SY(I)-SY(J) 2016 DZ=SZ(I)-SZ(J) 2017C call PBC (DX,DY,DZ) 2018 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 2019 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 2020 if(LHEX)then 2021C hexagonal periodic cell) 2022 XY=BOXYC*DX 2023 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 2024 DY=DY-BOYL 2025 DX=DX-HBOXL 2026 XY=BOXYC*DX 2027 end if 2028 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 2029 DY=DY+BOYL 2030 DX=DX+HBOXL 2031 XY=BOXYC*DX 2032 end if 2033 if(DY.gt.BOXY3+XY)then 2034 DY=DY-BOYL 2035 DX=DX+HBOXL 2036 XY=BOXYC*DX 2037 end if 2038 if(DY.lt.-BOXY3+XY)then 2039 DY=DY+BOYL 2040 DX=DX-HBOXL 2041 end if 2042 else 2043C rectangular cell 2044 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 2045 end if 2046 if(LOCT)then 2047C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 2048 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 2049 DX=DX-sign(CORSS,DX) 2050 DY=DY-sign(CORSS,DY) 2051 DZ=DZ-sign(CORSS,DZ) 2052 end if 2053 RR2=DX**2+DY**2+DZ**2 2054 if(LMNBL)then 2055 NMI=NSITS(ITYP) 2056 NMJ=NSITS(ITYPE(J)) 2057 if(NMI.gt.MMOL.and.NMJ.gt.MMOL)then 2058 RRM=RR2 2059 go to 130 2060 else if(NMI.gt.MMOL.and.NMJ.le.MMOL)then 2061 DX=SX(I)-X(JMOL) 2062 DY=SY(I)-Y(JMOL) 2063 DZ=SZ(I)-Z(JMOL) 2064 else if(NMJ.gt.MMOL)then 2065 DX=X(IMOL)-SX(J) 2066 DY=Y(IMOL)-SY(J) 2067 DZ=Z(IMOL)-SZ(J) 2068 else 2069 DX=X(IMOL)-X(JMOL) 2070 DY=Y(IMOL)-Y(JMOL) 2071 DZ=Z(IMOL)-Z(JMOL) 2072 end if 2073* call PBC (DX,DY,DZ) 2074 if(abs(DX) > HBOXL) DX = DX - sign(BOXL,DX) 2075 if(abs(DZ) > HBOZL) DZ = DZ - sign(BOZL,DZ) 2076 if(LHEX)then 2077C hexagonal periodic cell) 2078 XY=BOXYC*DX 2079 if(DY.gt.BOXY3-XY.and.(DX.gt.0..or.DY.gt.2.*BOXY3-XY))then 2080 DY=DY-BOYL 2081 DX=DX-HBOXL 2082 XY=BOXYC*DX 2083 end if 2084 if(DY.lt.-BOXY3-XY.and.(DX.lt.0..or.DY.lt.-2.*BOXY3-XY))then 2085 DY=DY+BOYL 2086 DX=DX+HBOXL 2087 XY=BOXYC*DX 2088 end if 2089 if(DY.gt.BOXY3+XY)then 2090 DY=DY-BOYL 2091 DX=DX+HBOXL 2092 XY=BOXYC*DX 2093 end if 2094 if(DY.lt.-BOXY3+XY)then 2095 DY=DY+BOYL 2096 DX=DX-HBOXL 2097 end if 2098 else 2099C rectangular cell 2100 if(abs(DY) > HBOYL) DY = DY - sign(BOYL,DY) 2101 end if 2102 if(LOCT)then 2103C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 2104 CORSS=HBOXL*int((4./3.)*(abs(DX)+abs(DY)+abs(DZ))/BOXL) 2105 DX=DX-sign(CORSS,DX) 2106 DY=DY-sign(CORSS,DY) 2107 DZ=DZ-sign(CORSS,DZ) 2108 end if 2109C square distance betwen COM of small molecules 2110 RRM=DX**2+DY**2+DZ**2 2111 else 2112 RRM=RR2 2113 end if 2114 130 if(IMOL.eq.JMOL)then 2115 if(RR2.le.RMX2)then 2116C Check that the atoms are non-bound 2117 do K=1,NNBB(I) 2118C 1-3 bound 2119 if(INBB(I,K).eq.J)go to 135 2120 if(INBB(I,K).eq.-J)then 2121C 1-4 bound 2122 MBSH=MBSH+1 2123 if(MBSH.gt.NBSMAX)write(*,*) 2124 + '!!! list of closest neigbours overfilled. atoms ',I,J 2125 NBS1(MBSH)=-I 2126 NBS2(MBSH)=J 2127 go to 135 2128 end if 2129 end do 2130 end if ! RR2.le 2131 end if ! IMOL.eq.JMOL 2132 if(RRM.le.SHORT2)then 2133 MBSH=MBSH+1 2134 if(MBSH.gt.NBSMAX)then 2135 write(*,*) 2136 +'!!! list of closest neigbours overfilled, atoms ',I,J,TASKID 2137 call FINAL 2138 end if 2139 NBS1(MBSH)=J 2140 NBS2(MBSH)=I 2141 else if(RRM.le.RSSQ)then 2142 MBLN=MBLN+1 2143 if(MBLN.gt.NBLMAX)then 2144 if(IOVER.eq.0)write(*,*) 2145 + '!!! list of far neigbours overfilled',TASKID 2146 IOVER=1 2147 end if 2148 if(IOVER.eq.0)then 2149 NBL1(MBLN)=J 2150 NBL2(MBLN)=I 2151 end if 2152 end if 2153C RDF calculation 2154 135 continue 2155*X IF(LGRC.and.RR2.le.RDFCUT2) THEN !=====> LGR 2156*X ISB = NSITE(J) 2157*X JSB = NSITE(I) 2158*X if(NGRI(ISB).le.NGRI(JSB))then 2159*X do JS=1,NGRI(ISB) 2160*X if(IGRI(ISB,JS).eq.JSB)then 2161*X IP = MGRI(ISB,JS) 2162*X go to 137 2163*X end if 2164*X end do 2165*X else 2166*X do JS=1,NGRI(JSB) 2167*X if(IGRI(JSB,JS).eq.ISB)then 2168*X IP = MGRI(JSB,JS) 2169*X go to 137 2170*X end if 2171*X end do 2172*X end if 2173*X go to 138 2174*X 137 R1 = sqrt(RR2) 2175*X GG = R1*DRDFI 2176*X MM = IDINT(GG)+1 2177*X if(MM.le.MAX)IRDF(MM,IP) = IRDF(MM,IP)+1 2178*X end if ! of LGRC 2179 138 continue 2180 end do 2181 end do 2182*X if(LGRC)NRDF=NRDF+1 2183* 2184* 15.3 Check that list of neighbours is not overfilled 2185* 2186 if(IOVER.ne.0)then 2187 write(*,*)' Insrease NBLMAX in prcm.h for this job to ',MBLN 2188 I=MBLN/NTOT+1 2189 write(*,*)' (change to: NBLMAX =',I,'* NTOT)' 2190 call FINAL 2191 end if 2192 timev = timev+cputime(timev0) 2193 if(IPRINN.lt.7)return 2194 if(LPIMD)then 2195 write(*,*)' node ',JBEAD, 2196 +'total num.of interactions: ',MBSH,' and ',MBLN 2197 else 2198 write(*,*)' Task ',TASKID, 2199 +'total num.of interactions: ',MBSH,' and ',MBLN 2200 end if 2201 return 2202 end 2203 2204