1*===================================================================== 2* MDynaMix v.5.2 3* Part 3 4* 5* Simulation setup 6* ---------------- 7C 8C This file contains subroutines setting up simulation: 9C 10C 1. IUNITS - set up units, addresses, references 11C 2 COMB_RULE - combination rule for non-bonded LJ interactions 12C 2. BNDLST - create list of bound, nonbound and 13C 1-4 neighbouring atoms 14C 3. CUBIC - set up molecular COM on cubic lattice 15C 4. FCC - set up molecular COM on FCC lattice 16C 5. DNASOL - put some molecule(s) in a cylindrical hole and 17C distribute other around 18C 6. SPHOLE - put some molecule(s) in a spherical hole and 19C distribute other around 20C 7. MIXTURE - mix molecules randomly 21C 8. PULLINI - initialization of restraints 22C 9. MKPARLIST- create list of pairs for LJ interactions 23C 24*=============== IUNITM =============================================== 25* 26* 27 SUBROUTINE IUNITS 28 include "prcm.h" 29 character*128 TAKESTR,STR 30* 31* 1.1 Calculate box dimensions: 32* ---------------------------- 33* 1.1.1 Some constants 34 FACTM = AVSNO*1000.0 35 TOTMAS = TOTMAS/FACTM ! in kg 36 FNOP = DFLOAT(NOP) 37 AVOFNR = AVSNO/FNOP 38 IF(RHO.NE.0.D0) THEN 39 VOLM = AVOFNR*(TOTMAS/(RHO*1000.0)) 40* 1.1.2 Define box size from density 41 IF(BOXL*BOYL*BOZL.lt.1.D-20) THEN 42 if(TASKID.eq.MAST)then 43 write(*,*)' Define box size from the density' 44 write(*,*) 45 end if 46 VOL = VOLM*FNOP*1.d30/AVSNO 47 if(ICELL.eq.1)VOL=2.*VOL 48 if(ICELL.eq.2)VOL=VOL/cos30 49 CL = VOL**(1.0/3.0) 50 BOXL=CL 51 BOYL=CL 52 BOZL=CL 53C Truncated octahedron geometry 54 if(ICELL.eq.1)VOL=0.5*VOL 55C Hexagonal cell geometry 56 if(ICELL.eq.2)then 57 BOYL=BOXL*cos30 58 VOL=BOXL*BOYL*BOZL 59 end if 60* 1.1.3 Box size difened from the input file 61 else 62 CL=BOXL 63 ENDIF 64 ELSE 65* 1.1.4 Case of vacuum simulation 66 BOXL = 1000.D0 67 BOYL = 1000.D0 68 BOZL = 1000.D0 69C No treatment of long-range electrostatics 70 ALPHA = 0.d0 71 if(TASKID.eq.MAST)PRINT "(/,'*** VACUUM SIMULATION ***',/)" 72 END IF 73C remember size/density from input 74C ( they will be owerwritten by the restart file) 75 BOXLT=BOXL 76 BOYLT=BOYL 77 if(LHEX)BOYLT=BOXLT*cos30 78 BOZLT=BOZL 79 RHON=RHO 80 if(ICELL.eq.0.and.TASKID.eq.MAST) 81 +write(*,*)' Rectangular geometry is used' 82 if(ICELL.eq.1.and.TASKID.eq.MAST) 83 +write(*,*)' Truncated octahedron geometry is used' 84 if(ICELL.eq.2.and.TASKID.eq.MAST)write(*,*) 85 +' Hexagonal periodic cell is used' 86* 1.1.5 Calculate size-dependent variables 87 call RECLEN(ODIN,ODIN,ODIN) 88 VOLM = 1.d-30*VOL*AVOFNR 89 DT = DT*1.d-15 ! in sec 90 TTREK=TTREK*1.d-15 91* 92* 1.2 Define different constants 93* ---------------------------- 94* 1.2.1 Internal units: 95 UNITL = 1.d-10 ! m 96 UNITM = TOTMAS ! mass of the system in kg 97 UNITT = DT ! time step 98 UNITE = UNITM*(UNITL/UNITT)**2 99 UNITP = UNITE*1.e-5/UNITL**3 100* 1.2.2 Dimensionless quantities and auxilarly constants 101 TSTEP = UNITT/SQRT((UNITM/UNITE)*UNITL**2) 102 HTSTEP = 0.5*TSTEP 103 TSTEB = TSTEP*FTIMES 104 HTSTEB = 0.5*TSTEB 105 BETA = UNITE/(BOLTZ*TRTEMP) 106 ENERF = AVSNO*UNITE 107 TOKJ = 0.001*ENERF 108 TEMPF = 2.0/(3.0*AVSNO*BOLTZ) 109 COULF = (ELECHG/EPS0)*(ELECHG/UNITE)/(4.D0*PI*UNITL) 110 PERMOL = ENERF*1.D-3/FNOP 111 SHORT2 = SHORT**2 112 SHORTA = SHORT-1.d-11/UNITL ! -0.1A 113 RLR = 1.d0/(BETA*RLR**2) 114 LSCFL = .false. 115 RAEX = 0. 116* 1.2.3 Some output 117 if(IPRINT.ge.7)then 118 PRINT "(/'*** UNIT MASS ',D13.6,' kg ' )",UNITM 119 PRINT "( '*** UNIT TIME ',D13.6,' S ' )",UNITT 120 PRINT "( '*** UNIT ENERGY ',D13.6,' J ' )",UNITE 121 PRINT "( '*** UNIT LENGTH ',D13.6,' m ' )",UNITL 122 PRINT"('*** UNIT PRESSURE ',D13.6,' 10**5 N/m**2 '/)",UNITP 123 PRINT "( '*** ENERGY FACTOR ',D13.6,' J/mol ' )",ENERF 124 PRINT "( '*** TEMP FACTOR ',D13.6,' K/J ' )",TEMPF 125 PRINT "( '*** COULOMB FACTOR ',D13.6,' '/)",COULF 126 PRINT "( '*** 1/kT (beta) ',D13.6,' u.e. ' )",BETA 127 PRINT "( '*** LONG TIME STEP ',F13.6,' .I.U. ' )",TSTEP 128 if(.not.LSHEJK)PRINT 129 +"( '*** SHORT TIME STEP ',F13.6,' .I.U. ' )",TSTEB 130 end if 131* 132* 1.2.4 Calculate energy/temperature convertion factors 133* --------------------------------------------------- 134C These factors are calculated for each molecule type taking into 135C account possible constraints. 136 CONVET = 3.*ENERF*TEMPF/FNST 137 do ITYP=1,NTYPES 138 FNSS3 = 3.*DFLOAT(NSITS(ITYP)) 139 FNSS1 = NSTOT/(NSTOT-1.) 140 TFACT(ITYP) = ENERF*TEMPF*FNSS1/(NSITS(ITYP)*NSPEC(ITYP)) 141 if(LSHEJK.and.ISHEJK(ITYP).ne.0)then 142 TFACT(ITYP)=TFACT(ITYP)*FNSS3/(FNSS3-NRCON(ITYP)) 143 end if 144 NRTSC(ITYP) = 0 ! will report num. of temperature scalings 145 end do 146* 147* 1.2.5 dimensionless mass parameters: 148 DO IS = 1,NSITES 149 MASSD (IS) = MASS (IS)/(NBEADS*FACTM*UNITM) ! fraction of the total mass 150 END DO! OF I 151* 152* 1.3 Addresses and references 153* ---------------------------- 154* 1.3.1 Addresses arrays 155C See comment to p.4 in main.f 156C The addresses provide atom number corresponding 157C to each new molecule, site ,type 158 IATOM = 0 159 IMOL = 0 160 IST = 0 161 IADDR (1) = 0 162 ISADR (1) = 0 163 ISADDR(1) = 0 164 TOTCH = 0. 165 NTYP1 = NTYPES + 1 166 DO ITYP = 1,NTYPES 167 IADDR (ITYP+1) = IADDR (ITYP)+NSPEC(ITYP) ! molecules 168 ISADR (ITYP+1) = ISADR (ITYP)+NSITS(ITYP) ! sites 169 ISADDR(ITYP+1) = ISADDR(ITYP)+NSITS(ITYP)*NSPEC(ITYP) ! atoms 170* 171* 1.3.2 Setup reference arrays 172* -------------------------- 173C Now, we can create list of references 174C 175C Reference arrays 176C NNUM atom -> molec 177C NSITE atom -> site 178C ITYPE atom -> type 179C ITS site -> type 180C ITM molec -> type 181 if(TASKID.eq.MAST) 182 + write(*,'(a6,i4,a3,a16,3x,a5,i5,4x,a6,i4)') 183 + ' type ',ITYP,' - ',NAMOL(ITYP)(1:16),' Nmol',NSPEC(ITYP), 184 + 'Nsites',NSITS(ITYP) 185 do IS=1,NSPEC(ITYP) 186 IMOL=IMOL+1 187 do I=1,NSITS(ITYP) 188 ISITE=IST+I 189 IATOM=IATOM+1 190 NNUM(IATOM)=IMOL ! atom -> molec 191 NSITE(IATOM)=ISITE ! atom -> site 192 ITYPE(IATOM)=ITYP ! atom -> type 193 ITS(ISITE)=ITYP ! site -> type 194 if(IPRINT.ge.10)write(*,'(4I7,3x,2a4,a6)') 195 + IATOM,IMOL,ISITE,ITYP,NM(ISITE),' on ',NAME(ITYP) 196 TOTCH = TOTCH+CHARGE(ISITE) 197 end do 198 ITM(IMOL)=ITYP ! molec -> type 199 end do 200 IST=IST+NSITS(ITYP) 201 END DO! OF ITYP 202* 203* 1.4 Report some quantities 204* --------------------------- 205 if(TASKID.eq.MAST)then 206 write(*,*) 207 write(*,'(a,f8.3)')'*** Total charge Q = ',TOTCH 208 end if 209* 210 if(IPRINT.ge.10)then 211 write(*,*) 212 PRINT "(25X,'*** ADDRESSES: ***')" 213 DO ITYP = 1,NTYP1 214 PRINT "(10X,4(1X,I8))", 215 + ITYP,ISADR(ITYP),IADDR(ITYP),ISADDR(ITYP) 216 END DO! OF ITYP 217 end if 218* 219* 1.5 Setup values of different working arrays 220* ------------------------------------------- 221 N = 0 222 MT = 0 223* 1.5.1 Inverse masses, charges 224 DO ITYP = 1,NTYPES 225 QS = 0.D0 226 NSB = ISADR (ITYP)+1 227 NSE = ISADR (ITYP+1) 228 DO I = 1,NSPEC(ITYP) 229 DO IS = NSB,NSE 230 N = N+1 231 MASSDI(N) = 1.0/MASSD(IS) 232 Q (N) = CHARGE(IS) 233 QS = QS+DABS(Q(N)) 234 END DO! OF I 235 END DO! OF IS 236 QSUM(ITYP) = QS 237* 1.5.2 Type pair indeces; energy components 238 DO JTYP = ITYP,NTYPES 239 MT = MT+1 240 MDX(ITYP,JTYP) = MT 241 MDX(JTYP,ITYP) = MT 242 POTES(MT) = 0.D0 243 POTLJ(MT) = 0.D0 244 END DO! OF JTYP 245* Drag g-force - conversion from 10^12 m/s^2 (1 m/s in 1 ps) to internal units 246 if(LDRAGT(ITYP))then 247 write(*,*)' g-force on type ',ITYP, 248 + ' : ',FDRAG(ITYP),' 10^12 m/s^2' 249 end if 250 FDRAG(ITYP)=FDRAG(ITYP)*(1.d12*UNITT**2/UNITL) 251 END DO! OF ITYP 252 MOLINT = MT+1 253 POTES(MOLINT)=0. 254 POTLJ(MOLINT)=0. 255C This counts energy absorbed from the thermostat 256 EABS = 0. 257 EABSS = 0. 258* 1.5.3 External field parameters 259C Frequence and amplitude of the external field in int.units 260 if(IEXT.eq.1)then 261 write(*,'(a,e12.5,a,e12.5,a)')' External filed: ',EXAMPL, 262 +' V/cm; Frequence ',EXFREQ,' hz' 263 EXAMPL=EXAMPL*UNITL*1.d2*ELECHG/UNITE ! V/cm -> i.u. 264 write(*,*)' In internal units: ',EXAMPL 265 end if 266* 267 if(IPRINT.ge.9) 268 +PRINT "(/' *** ',I2,' MOLECULAR INTERACTIONS ***'/)",MT 269* 270 DO ITYP = 1,NTYPES 271 DO JTYP = ITYP,NTYPES 272 MT = MDX(ITYP,JTYP) 273 if(IPRINT.ge.9) 274 +PRINT "(' NR: ',I2,5X,A6,' <--> ',A6)",MT,NAME(ITYP),NAME(JTYP) 275 END DO! OF JTYP 276 END DO! OF ITYP 277 LMOL1=.false. 278 LMOL2=.false. 279* 280* 1.6 Creating bonds, angles, torsion lists - for each node 281* --------------------------------------------------------- 282* 1.6.1 Loop over molecule types 283 NRBT=0 284 NRAT=0 285 NRRT=0 286 NRIT=0 287 ICNB=0 288 ICNA=0 289 ICNT=0 290 if(LBONDL)then 291 open(unit=26,file="bond.list",status='unknown') 292 write(26,'(a,a)')'# Bond list for job ',FNAME 293 end if 294 do ITYP=1,NTYPES 295 NSS=NSPEC(ITYP) 296* 1.6.2 Create bond list 297 if(IPOT(ITYP).eq.0)then ! otherwise special subroutine used (STRBND) 298 NBS=NRB(ITYP) 299 JSFB=IADB(ITYP)-1 300 NAS=NRA(ITYP) 301 JSFA=IADA(ITYP)-1 302 NTS=NRT(ITYP) 303 JSFT=IADT(ITYP)-1 304 do I=1,NSS ! over molecules 305 IML=ISADDR(ITYP)+(I-1)*NSITS(ITYP) 306 do J=1,NBS ! over bonds 307 IS=JSFB+J 308 IBD=ICNB+(I-1)*NBS+J 309 if(mod(IBD,NUMTASK).eq.TASKID)then 310 NRBT=NRBT+1 311 if(NRBT.gt.NBO)stop "increase NBO" 312 IBI(NRBT)=IB(IS)+IML 313 JBJ(NRBT)=JB(IS)+IML 314 IBUK(NRBT)=IS 315 end if 316C write bond list to file "bonds.list" if specified 317 if(LBONDL)then 318 IS1=IB(IS)+ISADR(ITYP) 319 IS2=JB(IS)+ISADR(ITYP) 320 write(26,'(2I6,6x,A4,A3,2A4,A6,2x,3I5,a)') 321 + IB(IS)+IML,JB(IS)+IML,NM(IS1),' - ',NM(IS2),' at ',NAME(ITYP),I 322 end if 323 end do 324* 1.6.3 Create covalent angle list 325 do J=1,NAS 326 IS=JSFA+J 327 IBD=ICNA+(I-1)*NAS+J 328 if(mod(IBD,NUMTASK).eq.TASKID)then 329 NRAT=NRAT+1 330 if(NRAT.gt.NAO)stop "increase NAO" 331 IAI(NRAT)=IA(IS)+IML 332 JAJ(NRAT)=JA(IS)+IML 333 KAK(NRAT)=KA(IS)+IML 334 IAUK(NRAT)=IS 335 end if 336 end do 337* 1.6.4 Create dihedral angle list 338 do J=1,NTS 339 IS=JSFT+J 340 IBD=ICNT+(I-1)*NTS+J 341 if(mod(IBD,NUMTASK).eq.TASKID)then 342 NRRT=NRRT+1 343 if(NRRT.gt.NTO)stop "increase NTO" 344 ITI(NRRT)=IT(IS)+IML 345 JTJ(NRRT)=JT(IS)+IML 346 KTK(NRRT)=KT(IS)+IML 347 LTL(NRRT)=LT(IS)+IML 348 ITUK(NRRT)=IS 349 end if 350 end do 351 end do 352 ICNB=ICNB+NSS*NBS 353 ICNA=ICNA+NSS*NAS 354 ICNT=ICNT+NSS*NTS 355 else 356* 1.6.5 Generate list of bonds for special potential model 357C (only for output; not used in MD) 358 if(LBONDL)then 359 NBS=NRB(ITYP) 360 JSFB=IADB(ITYP)-1 361 do I=1,NSS ! over molecules 362 IML=ISADDR(ITYP)+(I-1)*NSITS(ITYP) 363 do J=1,NBS ! over bonds 364 IS=JSFB+J 365 IS1=IB(IS)+ISADR(ITYP) 366 IS2=JB(IS)+ISADR(ITYP) 367 write(26,'(2I6,6x,A4,A3,2A4,A6,2x,3I5,a)') 368 + IB(IS)+IML,JB(IS)+IML,NM(IS1),' - ',NM(IS2),' at ',NAME(ITYP),I 369 end do 370 end do 371 end if 372 end if 373 end do 374C write(*,*)' node ',TASKID,' angles ',NRAT 375C do I=1,NRAT 376C write(*,*)TASKID,I,IAI(I),JAJ(I),KAK(I),IAUK(I) 377C end do 378C write(*,*)' node ',TASKID,' torsions ',NRRT 379C do I=1,NRRT 380C write(*,*)TASKID,I,ITI(I),JTJ(I),KTK(I),LTL(I),ITUK(I) 381C end do 382* 1.6.5 Calculate lengths of transfer arrays 383 LENMB=8*NRBB 384 LENMA=8*NRAA 385 LENMT=8*NRTT 386 LENMP=8*NSTOT 387 LENIT=8*MOLINT 388 LENTT=8*NTYPES 389C In unused.f : 390C Additional list of bonds 391C call ADDBONDS 392* 393* 1.7 CONVERSION TO INTERNAL UNITS: 394* ---------------------------------- 395* 1.7.1 Energy factors 396 FACTOR = 1.D3/AVSNO/UNITE 397* EFACT = FKCALJ*FACTOR ! kcal 398 EFACT = FACTOR ! kJ 399* 1.7.2 Bonds force cnstants 400C Harmonic bonds: 401 DO I = 1,NB 402 FB(I) = FB(I)*EFACT 403 END DO! OF I 404C Morse potentials: 405 DO I = 1,NB 406 DB(I) = DB(I)*EFACT 407 END DO! OF I 408* 409* 1.7.3 Angle bending: 410* 411 DO I = 1,NA 412 FA(I) = FA(I)*EFACT 413 RA(I) = RA(I)/TODGR 414 END DO! OF I 415* 416* 1.7.4 Torsional angles 417* 418 DO I = 1,NT 419 RT (I) = RT (I)/TODGR 420 FT (I) = FT (I)*EFACT 421 FT1(I) = FT1(I)*EFACT 422 FT2(I) = FT2(I)*EFACT 423 FT3(I) = FT3(I)*EFACT 424 FT4(I) = FT4(I)*EFACT 425 FT5(I) = FT5(I)*EFACT 426 END DO! OF I 427* 428* 1.8 Setting up non-bonded LJ interactions 429* ----------------------------------------- 430* 1.8.1 Assigment of non-bonded interaction types 431* 1.8.1.1 432* If we use special list of LJ cross-pairs, interaction types = sites 433 if(LPAIRS.or.NNADP.gt.0)then 434 do I=1,NSITES 435 SIGT(I)=SIGMA(I) 436 EPST(I)=EPSIL(I) 437 SIGAT(I)=SIGAD(I) 438 EPSAT(I)=EPAD(I) 439 INBT(I)=I 440 end do 441 NNBT = NSITES 442* 1.8.1.2 sites with same LJ parameters and charges are assigned the same interaction type 443 else ! .not.LPAIRS 444 NNBT=1 445 SIGT(1)=SIGMA(1) 446 EPST(1)=EPSIL(1) 447 do J=1,NNAD 448 if(ILJ(J).eq.1)then 449 SIGAT(1)=SIGAD(J) 450 EPSAT(1)=EPAD(J) 451 go to 181 452 end if 453 end do 454 SIGAT(1)=SIGMA(1) 455 EPSAT(1)=EPSIL(1) 456 181 continue 457 INBT(1)=1 458 do IS=2,NSITES 459 ITYP = ITS(IS) 460 do JS=1,NNBT 461 if(dabs(SIGMA(IS)-SIGT(JS)).lt.1.d-12.and. 462 + dabs(EPSIL(IS)-EPST(JS)).lt.1.d-12)then 463* found same eps and sigma as interaction type JS 464C check for EE sites (expanded ensemble ) 465 if(IEE(ITYP).eq.0)then 466 do J=1,IS-1 467 JTYP=ITS(J) 468 JSS=INBT(J) 469 if(JSS.eq.JS.and.IEE(JTYP).ne.0)go to 183 470 end do 471 end if 472 if(IEE(ITYP).eq.1)then 473 do J=1,IS-1 474 JTYP=ITS(J) 475 JSS=INBT(J) 476 if(JSS.eq.JS.and.IEE(JTYP).eq.0)go to 183 477 end do 478 end if 479* check for 1-4 480 do J=1,NNAD 481 if(ILJ(J).eq.IS)go to 182 482 end do 483 if(dabs(SIGMA(IS)-SIGAT(JS)).lt.1.d-12.and. 484 + dabs(EPSIL(IS)-EPSAT(JS)).lt.1.d-12)then 485 INBT(IS)=JS 486 go to 186 ! yes 487 else 488 go to 183 ! no 489 end if 490 182 continue 491 if(dabs(SIGAD(J)-SIGAT(JS)).lt.1.d-12.and. 492 + dabs(EPAD(J)-EPSAT(JS)).lt.1.d-12)then 493* this site (IS) has a nonbonded type JS 494 INBT(IS)=JS 495 go to 186 496 end if 497 end if 498 end do 499* no previous site found - set up new non-bonded type 500 183 NNBT=NNBT+1 501 if(NNBT.gt.NNBTM)STOP "STOP: Increase NNBTM in dimpar.h !!!" 502 INBT(IS)=NNBT 503 SIGT(NNBT)=SIGMA(IS) 504 EPST(NNBT)=EPSIL(IS) 505* check for non-standard 1-4 LJ parameters 506 do J=1,NNAD 507 if(ILJ(J).eq.IS)then 508 SIGAT(NNBT)=SIGAD(J) 509 EPSAT(NNBT)=EPAD(J) 510 go to 186 511 end if 512 end do 513 SIGAT(NNBT)=SIGT(NNBT) 514 EPSAT(NNBT)=EPST(NNBT) 515 186 continue 516 end do 517 end if ! LPAIRS 518 if(TASKID.eq.MAST) 519 + write(*,*)NNBT,' different non-bonded types found' 520 if(IPRINT.ge.7.and.TASKID.eq.MAST)then 521 write(*,*)' Non-bonded types:' 522 write(*,*) 523 +' site typ typ_num iee eps sig eps14 sig14' 524 do I=1,NSITES 525 ITYP=ITS(I) 526 J=INBT(I) 527 write(*,'(4i6,4f12.4)') 528 + I,ITYP,J,IEE(ITYP),SIGT(J),EPST(J),SIGAT(J),EPSAT(J) 529 end do 530 end if 531* 1.8.2 Setting up LJ cross-interactions 532C array EPST to internal energy units 533 do IS=1,NNBT 534 EPST(IS)=EFACT*EPST(IS) 535 EPSAT(IS)=EFACT*EPSAT(IS) 536 end do 537 do IS=1,NNBT 538 do JS=1,NNBT 539* non-bonded LJ 540 EPS1=EPST(IS) 541 EPS2=EPST(JS) 542 call COMB_RULES(SIGT(IS),SIGT(JS),EPS1,EPS2,A6,B12,ICR) 543 A6LJ(IS,JS)=A6 544 B12LJ(IS,JS)=B12 545* "0" might be reserved for intramolecular; without "0" may be subject to change 546* in expanded ensemble 547 A6LJ0(IS,JS)=A6 548 B12LJ0(IS,JS)=B12 549* 1-4 LJ 550 EPS1=EPSAT(IS) 551 EPS2=EPSAT(JS) 552 call COMB_RULES(SIGAT(IS),SIGAT(JS),EPS1,EPS2,A6,B12,ICR) 553 A6_14(IS,JS)=A6 554 B12_14(IS,JS)=B12 555 end do 556 end do 557* 1.8.3 assigning special 1-4 pairs from .mmol files 558 do IP=1,NNADP 559 IS = ILJP(IP) ! site = int.type 560 JS = JLJP(IP) 561 A6 = -4.d0*EFACT*EPAP(IP)*SIGAP(IP)**6 562 B12 = 4.d0*EFACT*EPAP(IP)*SIGAP(IP)**12 563 A6_14(IS,JS) = A6 564 B12_14(IS,JS) = B12 565 A6_14(JS,IS) = A6 566 B12_14(JS,IS) = B12 567 end do 568* 1.8.4 reading special cross-term pairs 569 if(LPAIRS)then 570 if(LPARFIL)then 571 call MKPARLIST 572 else 573 open(unit=23,file=fpairs,status='old',err=198) 574 end if 575 if(TASKID.eq.MAST)write(*,*) 576 +' Redefine LJ parameters for pairs: ' 577 191 IE=-20 578 STR=TAKESTR(23,IE) 579 if(IE.eq.-1)go to 197 580 read(STR,*,end=194,err=194)II,JJ,SIG,EPS,SIG14,EPS14 581 A6LJ(II,JJ)= -4.d0*EFACT*EPS*SIG**6 582 B12LJ(II,JJ)= 4.d0*EFACT*EPS*SIG**12 583 A6LJ0(II,JJ)= -4.d0*EFACT*EPS*SIG**6 584 B12LJ0(II,JJ)= 4.d0*EFACT*EPS*SIG**12 585 A6_14(II,JJ)= -4.d0*EFACT*EPS14*SIG14**6 586 B12_14(II,JJ)= 4.d0*EFACT*EPS14*SIG14**12 587 A6LJ(JJ,II)= -4.d0*EFACT*EPS*SIG**6 588 B12LJ(JJ,II)= 4.d0*EFACT*EPS*SIG**12 589 A6LJ0(JJ,II)= -4.d0*EFACT*EPS*SIG**6 590 B12LJ0(JJ,II)= 4.d0*EFACT*EPS*SIG**12 591 A6_14(JJ,II)= -4.d0*EFACT*EPS14*SIG14**6 592 B12_14(JJ,II)= 4.d0*EFACT*EPS14*SIG14**12 593 if(TASKID.eq.MAST)write(*,'(4a,2i5,a,4f10.4)') 594 + NM(II),'-',NM(JJ),' (',II,JJ,') :',SIG,EPS,SIG14,EPS14 595 go to 191 596 194 read(STR,*,end=197,err=197)II,JJ,SIG,EPS 597 A6LJ(II,JJ)= -4.d0*EFACT*EPS*SIG**6 598 B12LJ(II,JJ)= 4.d0*EFACT*EPS*SIG**12 599 A6LJ0(II,JJ)= -4.d0*EFACT*EPS*SIG**6 600 B12LJ0(II,JJ)= 4.d0*EFACT*EPS*SIG**12 601 A6LJ(JJ,II)= -4.d0*EFACT*EPS*SIG**6 602 B12LJ(JJ,II)= 4.d0*EFACT*EPS*SIG**12 603 A6LJ0(JJ,II)= -4.d0*EFACT*EPS*SIG**6 604 B12LJ0(JJ,II)= 4.d0*EFACT*EPS*SIG**12 605 if(TASKID.eq.MAST)write(*,'(4a,2i5,a,4f10.4)') 606 + NM(II),'-',NM(JJ),' (',II,JJ,') :',SIG,EPS 607 go to 191 608 197 close(23) 609 go to 199 610 198 if(TASKID.eq.MAST)then 611 write(*,*)' File with special LJ pairs not found' 612 write(*,*)' Usual combination rules will be used' 613 end if 614 199 continue 615 end if ! LPAIRS 616 if(IPRINT.ge.9.and.TASKID.eq.MAST)then 617 write(*,*)' A6s ' 618 do IS=1,NSITES 619 II=INBT(IS) 620 write(*,'(2i4,2000e12.4)') 621 +IS,II,(A6LJ(II,INBT(JS)),JS=1,NSITES) 622 end do 623 write(*,*)' B12s ' 624 do IS=1,NSITES 625 II=INBT(IS) 626 write(*,'(2i4,2000e12.4)') 627 +IS,II,(B12LJ(II,INBT(JS)),JS=1,NSITES) 628 end do 629 end if 630* 631* 1.9 Distribution of atoms over procesors 632* ---------------------------------------- 633C (Relevant only for parallel execution) 634 NAPC=NSTOT/NUMTASK+1 635 IF((NSTOT-((NUMTASK-1)*NAPC)).LT.0) THEN 636 NAPC = NAPC - 1 637 END IF 638 if(LSHEJK)then 639* 1.9.1 distribution of molecules - trying to give equal work for SHEJK 640 I=0 641 NABS(0)=0 642 NAB(0)=1 643 IWRK=0 644 AWRK=NRBND*1./NUMTASK 645 do ITYP=1,NTYPES 646 IBEG=IADDR(ITYP)+1 647 IEND=IADDR(ITYP+1) 648 NSS =NSITS(ITYP) 649 NWRK=NRB(ITYP) 650 if(ITYP.lt.NTYPES)NWRK1=NRB(ITYP+1) 651 do IMOL=IBEG,IEND 652 IWRK=IWRK+NWRK 653 DIFF=IWRK-AWRK*(I+1) 654 if(DIFF.ge.0)then 655 NAE(I)=ISADDR(ITYP)+(IMOL-IBEG)*NSS 656 if(DIFF.le.0.5*NWRK.or.NAE(I).eq.NABS(I))NAE(I)=NAE(I)+NSS 657 NAP(I)=NAE(I)-NABS(I) 658 NAP3(I)=3*NAP(I) 659 I=I+1 660 if(I.ge.NUMTASK)go to 65 661 NABS(I)=NAE(I-1) 662 NAB(I)=NABS(I)+1 663 end if 664 end do 665 end do 666 65 NAE(NUMTASK-1)=NSTOT 667 NAP(NUMTASK-1)=NSTOT-NABS(NUMTASK-1) 668 NAP3(NUMTASK-1)=3*NAP(NUMTASK-1) 669 else ! not.SHEJK 670* 1.9.2 Equilibrium distribution of atoms over processors 671 do I=0,NUMTASK-1 672 NABS(I)=NAPC*I 673 NAB(I)=NABS(I)+1 674 NAE(I)=NABS(I)+NAPC 675 NAP(I)=NAPC 676 NAP3(I)=3*NAPC 677 NABS3(I)=3*NABS(I) 678 end do 679 if(NAE(NUMTASK-1).gt.NTOT)stop 'increase NTOT for overhead' 680 NAE(NUMTASK-1)=NSTOT 681 NAP(NUMTASK-1)=NSTOT-NABS(NUMTASK-1) 682 NAP3(NUMTASK-1)=3*NAP(NUMTASK-1) 683 end if 684 do I=0,NUMTASK-1 685 NABS3(I)=3*NABS(I) 686 end do 687 if(IPRINT.ge.7)then 688 write(*,*)' node Nat from to NABS ' 689 do I=0,NUMTASK-1 690 write(*,'(5i7)')I,NAP(I),NAB(I),NAE(I),NABS(I) 691 end do 692 end if 693C Harmonic constant for PIMD 694 BEADFAC = 0.5d-3*NBEADS*(BOLTZ*TRTEMP/HPLANK)**2*UNITL**2 695 / /(UNITE*AVSNO) 696** write(*,*)JBEAD,' BEADFAC=',BEADFAC*UNITE/UNITL**2,' J/m^2' 697*** FBMASS=1.-1./FBMASS 698 CONVEQ = ENERF*TEMPF/NBEADS 699 RETURN 700 END 701* 702*================= COMB_RULES ==================================== 703* 704* 2. COMB_RULES: 705* Combination rules for LJ parameters 706* 707 subroutine COMB_RULES(S1,S2,E1,E2,A6,B12,ICR) 708 real*8 S1,S2,E1,E2,A6,B12 709* Geometrical rule for sigma 710 if(ICR.eq.1)then 711 A6 = -4.d0*sqrt(E1*E2)*(S1*S2)**3 712 B12 = 4.d0*sqrt(E1*E2)*(S1*S2)**6 713* Kong rules 714 else if(ICR.eq.2)then 715 A6 = -sqrt(E1*E2)*(S1*S2)**3 716 B12 =(E1*S1**12/2.d0**13)* 717 * (1+(E2*S2**12/(E1*S1**12))**(1./13.))**13 718* Fender_Halsey rules 719 else if(ICR.eq.3)then 720 A6 = -8.d0*(E1*E2)*(0.5*(S1+S2))**6/(E1+E2) 721 B12 = 8.d0*(E1*E2)*(0.5*(S1+S2))**12/(E1+E2) 722* Waldman-Hagler rules 723 else if(ICR.eq.4)then 724 SIP6 = (S1**6+S2**6)/2 725 EP = 2*sqrt(E1*E2)*(S1**3*S2**3/(S1**6+S2**6)) 726 A6 = -4.d0*EP*SIP6 727 B12 = 4.d0*EP*SIP6**2 728* Lorentz-Berthelot (default) 729 else 730 A6 = -4.d0*sqrt(E1*E2)*(0.5*(S1+S2))**6 731 B12 = 4.d0*sqrt(E1*E2)*(0.5*(S1+S2))**12 732 end if 733 return 734 end 735* 736*=============== BNDLST =========================================== 737* 738* 3. Prepare list of bonded atoms 739C 740 SUBROUTINE BNDLST 741* 742 include "prcm.h" 743 integer NRNOD(NTPS) 744* 745 do I=1,NSTOT 746 NNBB(I)=0 ! number of atoms bound to this 747 do IS=1,NBDMAX 748 INBB(I,IS)=0 ! which atoms are bound 749 end do 750 end do 751 IBEG=TASKID+1 752 NRN=NUMTASK 753 if(NRN.gt.NTYPES)NRN=NTYPES 754 do ITYP=1,NTYPES 755 NRNOD(ITYP)=mod(ITYP-1,NUMTASK) 756 end do 757 if(IBEG.gt.NTYPES)go to 101 758 do ITYP=IBEG,NTYPES,NUMTASK 759* 760 NSP = NSPEC(ITYP) ! num of mol of this type 761 NSS = NSITS(ITYP) ! num of sites on a mol of this type 762 NRBS = NRB(ITYP) ! num of bonds 763 NBBEG = IADB(ITYP) ! first bond 764 NBEND = NBBEG+NRBS-1 ! last bond 765 ISBEG = ISADR(ITYP)+1 ! first site 766 ISEND = ISADR(ITYP+1) ! last site 767 ISHF = ISADR(ITYP) ! shift of global site num. relatevely local 768* 769* First pass: 1-2 neigbours = bonds 770* 771 do I=NBBEG,NBEND 772 IBB=IB(I)+ISHF 773 JBB=JB(I)+ISHF 774 if(ID(I).ne.2)call BNDREG(IBB,JBB,1) 775 end do 776* 777* 1 - 3 neigbours: atom + its bonds 778* 779 do IS=ISBEG,ISEND 780 do J=NBBEG,NBEND 781 IBB=IB(J)+ISHF 782 JBB=JB(J)+ISHF 783 if(IS.eq.IBB.and.ID(J).ne.2)then 784 do K=J+1,NBEND 785 IBK=IB(K)+ISHF 786 JBK=JB(K)+ISHF 787 if(IS.eq.IBK.and.ID(K).ne.2)then 788 call BNDREG(JBB,JBK,1) 789 else if(IS.eq.JBK.and.ID(K).ne.2)then 790 call BNDREG(JBB,IBK,1) 791 end if 792 end do 793 else if(IS.eq.JBB.and.ID(J).ne.2)then 794 do K=J+1,NBEND 795 IBK=IB(K)+ISHF 796 JBK=JB(K)+ISHF 797 if(IS.eq.IBK.and.ID(K).ne.2)then 798 call BNDREG(IBB,JBK,1) 799 else if(IS.eq.JBK.and.ID(K).ne.2)then 800 call BNDREG(IBB,IBK,1) 801 end if 802 end do 803 end if 804 end do 805 end do 806* 807* 1 - 4 neighbours: bonds and its bonds 808* 1 - 4 neigbours are marked by minus in the list 809* 810 do I=NBBEG,NBEND 811 IS=IB(I) 812 JS=JB(I) 813 if(ID(I).ne.2)then 814 do J=NBBEG,NBEND 815 if(IS.eq.JB(J).and.ID(J).ne.2)then 816 do K=NBBEG,NBEND 817 if(JS.eq.IB(K).and.ID(K).ne.2)then 818 call BNDREG(IB(J)+ISHF,JB(K)+ISHF,-1) 819 else if(JS.eq.JB(K).and.ID(K).ne.2)then 820 call BNDREG(IB(J)+ISHF,IB(K)+ISHF,-1) 821 end if 822 end do 823 else if(IS.eq.IB(J).and.ID(J).ne.2)then 824 do K=NBBEG,NBEND 825 if(JS.eq.IB(K).and.ID(K).ne.2)then 826 call BNDREG(JB(J)+ISHF,JB(K)+ISHF,-1) 827 else if(JS.eq.JB(K).and.ID(K).ne.2)then 828 call BNDREG(JB(J)+ISHF,IB(K)+ISHF,-1) 829 end if 830 end do 831 end if 832 end do 833 end if 834 end do 835*---------------------------- 836 100 if(.not.LPIMD)then 837 write(*,*) 838 + ' node ',TASKID,' list of bound atoms for ',NAME(ITYP) 839 end if 840 end do ! of ITYP 841 101 call GATHER(NRNOD) 842CM 843 if(IPRINT.ge.8.and.TASKID.eq.MAST)then 844 write(*,*)' list of bound atoms: ' 845 do I=1,NSTOT 846 write(*,'(i5,i3,2x,40i5)')I,NNBB(I),(INBB(I,J),J=1,NNBB(I)) 847 end do 848 end if 849*---------------------------------------- 850* Check of special 1-4 pairs 851 do I=1,NNADP 852 IS = ILJP(I) 853 JS = JLJP(I) 854 ITYP = ITS(IS) 855 JTYP = ITS(JS) 856 if(ITYP.ne.JTYP)then 857 write(*,*)'!!! wrong types for special 1-4 pairs :' 858 write(*,*)' Sites: ',IS,JS,' Types: ',ITYP,JTYP 859 stop 860 end if 861 ISP = ISADDR(ITYP)+IS 862 JSP = ISADDR(ITYP)+JS 863 do J=1,NNBB(ISP) 864 if(JSP.eq.-INBB(ISP,J))go to 200 865 end do 866 write(*,'(a)') 867 + '!!! Warning: special 1-4 pairs are not 1-4 neigbours:' 868 write(*,'(a,i5,2(2x,a))') 869 + 'Type: ',ITYP,' Molecule: ',NAMOL(ITYP) 870 write(*,'(a,2i5,2(2x,a))') 871 + 'Sites: ',IS-ISADR(ITYP),JS-ISADR(ITYP),NM(IS),NM(JS) 872 200 continue 873 end do 874 RETURN 875 END 876* 877*============ BNDREG ========================================================== 878* 879* 3.2 registration of bound sites 880 subroutine BNDREG(IBB,JBB,ISIG) 881 include "prcm.h" 882 if(IBB.eq.JBB)return 883 ITYP=ITS(IBB) 884 NSP=NSITS(ITYP) 885 if(ISIG.ge.0)then 886 ISG=1 887 else 888 ISG=-1 889 end if 890 if(ITYP.ne.ITS(JBB))then 891 write(*,*)' Internal error: different types in BNDREG: ', 892 +IBB,JBB,':',ITYP,ITS(JBB) 893 stop 894 end if 895 ISHA = ISADDR(ITYP)-ISADR(ITYP) 896 do IM=1,NSPEC(ITYP) 897 IBS=IBB+ISHA+(IM-1)*NSP 898 JBS=JBB+ISHA+(IM-1)*NSP 899* check that this pair is not already registered 900 do I=1,NNBB(IBS) 901 if(iabs(INBB(IBS,I)).eq.JBS)go to 10 902 end do 903 NNBB(IBS)=NNBB(IBS)+1 904 INBB(IBS,NNBB(IBS))=ISG*JBS 905 if(NNBB(IBS).gt.NBDMAX)then 906 write(*,*) 907 +' list of bound atoms exceeded for atom ',IBS,' ',NM(IBS) 908 write(*,*)' Possibly, non-bonded interaction set off ???' 909 stop 910 end if 911 10 do I=1,NNBB(JBS) 912 if(iabs(INBB(JBS,I)).eq.IBS)go to 20 913 end do 914 915 NNBB(JBS)=NNBB(JBS)+1 916 INBB(JBS,NNBB(JBS))=ISG*IBS 917 if(NNBB(JBS).gt.NBDMAX)then 918 write(*,*) 919 +' list of bound atoms exceeded for atom ',JBS,' ',NM(JBS) 920 write(*,*)' Possibly, non-bonded interaction set off ???' 921 stop 922 end if 923 20 continue 924 end do 925 return 926 end 927* 928*================ CUBIC ================================================= 929* 930* 4. Set molecules on a cubic lattice 931* 932 SUBROUTINE CUBIC(X,Y,Z,BOXL,BOYL,BOZL,NOP,ICELL,IPRINT) 933 IMPLICIT real*8 (A-H,O-Z) 934* 935 DIMENSION X(*),Y(*),Z(*) 936* 937 if(NOP.eq.1)then 938 X(1)=0. 939 Y(1)=0. 940 Z(1)=0. 941 write(*,*)' Put molecule COM to 0 ' 942 return 943 end if 944 FNOP = DFLOAT(NOP) 945 if(ICELL.eq.1)then 946 FAC=2. 947 else 948 FAC=1. 949 end if 950 SIZE = (BOXL*BOYL*BOZL/FNOP/FAC)**(1.d0/3.d0) 951 NRX = BOXL/SIZE 952 NRY = BOYL/SIZE 953 NRZ = BOZL/SIZE 954 10 if(IPRINT.ge.7)write(*,*)NRX,NRY,NRZ 955 if(NRX*NRY*NRZ.lt.NOP)NRX=NRX+1 956 if(NRX*NRY*NRZ.lt.NOP)NRY=NRY+1 957 if(NRX*NRY*NRZ.lt.NOP)then 958 NRZ=NRZ+1 959 go to 10 960 end if 961 NRCUBE = NRX*NRY*NRZ 962 NEMPTY = NRCUBE-NOP 963 L = 0 964 DO 100 I = 1,NRX 965 DO 100 J = 1,NRY 966 DO 100 K = 1,NRZ 967 L = L+1 968 X(L) = (DFLOAT(I)-0.5)*BOXL/NRX-0.5*BOXL 969 Y(L) = (DFLOAT(J)-0.5)*BOYL/NRY-0.5*BOYL 970 Z(L) = (DFLOAT(K)-0.5)*BOZL/NRZ-0.5*BOZL 971 if(ICELL.eq.1.and.(abs(X(L))+abs(Y(L))+abs(Z(L))).gt.0.75*BOXL) 972 + L=L-1 973 if(L.ge.NOP)go to 110 974 100 CONTINUE 975 if(L.lt.NOP)then 976 NRX=NRX+1 977 go to 10 978 end if 979 110 NSKIP = 0 980 IF(NEMPTY.NE.0) NSKIP=NOP/NEMPTY+1 981 PRINT "(/1X,'*** CUBIC LATTICE: ',3I4,' SITES/EDGE -> ',I4/ 982 X' SITES TOTALLY ( ',I3,' VACANT - WITH INTERVALL OF ',I3,' )'/)" 983 X,NRX,NRY,NRZ,NRCUBE,NEMPTY,NSKIP 984* 985 DO I = 1,NOP 986 call PBC(X(I),Y(I),Z(I)) 987 end do 988* 989 call MIXTURE(X,Y,Z,NOP) 990* 991 RETURN 992 END 993* 994*=============== FCC =================================================== 995* 996* 5. Set molecules on the FCC lattice 997* ----------------------------------- 998 SUBROUTINE FCC(X,Y,Z,BOXL,BOYL,BOZL,NSP,ICELL) 999 IMPLICIT real*8 (A-H,O-Z) 1000 DIMENSION X(NSP),Y(NSP),Z(NSP) 1001 DIMENSION XB(4),YB(4),ZB(4) 1002 DATA XB/0.5,-.5,0.5,-.5/,YB/0.5,-.5,-.5,0.5/,ZB/0.5,0.5,-.5,-.5/ 1003C SET UP STARTING FCC LATTICE 1004 if(NSP.eq.1)then 1005 X(1)=0. 1006 Y(1)=0. 1007 Z(1)=0. 1008 write(*,*)' Put molecule COM to 0 ' 1009 return 1010 end if 1011 HBOXL = 0.5*BOXL 1012 VOL=BOXL*BOYL*BOZL 1013 if(ICELL.eq.1)then 1014 FAC=2. 1015 else 1016 FAC=1. 1017 end if 1018 SCELL=(4.*VOL/NSP)**0.333333333333 1019 NCX=BOXL/SCELL+0.999 1020 NCY=BOYL/SCELL+0.999 1021 NCZ=BOZL/SCELL+0.999 1022 5 continue 1023 if(4*NCX*NCY*NCZ.lt.NSP)then 1024 NCX=NCX+1 1025 if(4*NCX*NCY*NCZ.lt.NSP)then 1026 NCY=NCY+1 1027 if(4*NCX*NCY*NCZ.lt.NSP)then 1028 NCZ=NCZ+1 1029 go to 5 1030 end if 1031 end if 1032 end if 1033 10 FACT = HBOXL/DFLOAT(NCX) 1034 SHIFX = BOXL/DFLOAT(NCX) 1035 SHIFY = BOYL/DFLOAT(NCY) 1036 SHIFZ = BOZL/DFLOAT(NCZ) 1037 BASE = HBOXL*(-1.D0+1./DFLOAT(NCX)) 1038 IL = 0 1039 DO 20 IB = 1,4 1040 ZS = BASE+ZB(IB)*FACT 1041 DO 19 IZ = 1,NCZ 1042 YS = BASE+YB(IB)*FACT 1043 DO 18 IY = 1,NCY 1044 XS = BASE+XB(IB)*FACT 1045 DO 17 IX = 1,NCX 1046 IL = IL+1 1047 X(IL) = XS 1048 Y(IL) = YS 1049 Z(IL) = ZS 1050 if(ICELL.eq.1.and.(abs(X(IL))+abs(Y(IL))+abs(Z(IL))).gt. 1051 + 0.75*BOXL)IL=IL-1 1052 if(IL.ge.NSP)go to 24 1053 XS = XS+SHIFX 1054 17 CONTINUE 1055 YS = YS+SHIFY 1056 18 CONTINUE 1057 ZS = ZS+SHIFZ 1058 19 CONTINUE 1059 20 CONTINUE 1060 if(IL.lt.NSP)then 1061 IL=0 1062 NCY=NCY+1 1063 go to 10 1064 end if 1065* 1066 24 do I=1,NSP 1067 call PBC(X(I),Y(I),Z(I)) 1068 end do 1069* 1070 call MIXTURE(X,Y,Z,NSP) 1071* 1072 RETURN 1073 END 1074* 1075*=================== DNASOL ==================================== 1076* 1077* 6. Put one large molecule in a cylindrical hole along Z axis and 1078* distribute other molecules around it 1079* 1080 subroutine DNASOL 1081 include "prcm.h" 1082 dimension XRS(NPART),YRS(NPART),ZRS(NPART) 1083 NOPSM=0 1084 do ITYP=1,NTYPES 1085 if(IINIT(ITYP).ne.1)NOPSM=NOPSM+NSPEC(ITYP) 1086 end do 1087 VOLAV=VOL-BOZL*PI*RDNA**2 1088 if(TASKID.eq.MAST) 1089 +write(*,*)' Distribute ',NOPSM,' molecules in ',VOLAV,' A**3' 1090 SHAG=(VOLAV/NOPSM)**0.333333333 1091 IAT=0 1092 NSX=BOXL/SHAG+0.7 1093 NSY=BOYL/SHAG+0.7 1094 NSZ=BOZL/SHAG+0.5 1095 10 NCOUNT=0 1096 do IX=1,NSX 1097 XX=-HBOXL+BOXL*(dfloat(IX)-0.6)/NSX 1098 do IY=1,NSY 1099 YY=-HBOYL+BOYL*(dfloat(IY)-0.45)/NSY 1100 RR=sqrt(XX**2+YY**2) 1101 if(RR.gt.RDNA)then 1102 do IZ=1,NSZ 1103 NCOUNT=NCOUNT+1 1104 XRS(NCOUNT)=XX 1105 YRS(NCOUNT)=YY 1106 ZRS(NCOUNT)=-HBOZL+BOZL*(dfloat(IZ)-0.45)/NSZ 1107 if(NCOUNT.ge.NOPSM)go to 20 1108 end do 1109 end if 1110 end do 1111 end do 1112 if(NCOUNT.le.NOPSM)then 1113 IAT=IAT+1 1114 if(mod(IAT,3).eq.0)then 1115 NSX=NSX+1 1116 if(TASKID.eq.MAST)write(*,*)' increase NSX to ',NSX 1117 else if(mod(IAT,3).eq.1)then 1118 NSY=NSY+1 1119 if(TASKID.eq.MAST)write(*,*)' increase NSY to ',NSY 1120 else 1121 NSZ=NSZ+1 1122 if(TASKID.eq.MAST)write(*,*)' increase NSZ to ',NSZ 1123 end if 1124 go to 10 1125 end if 1126* peremeshivanie 1127 20 call MIXTURE(XRS,YRS,ZRS,NOPSM) 1128 J=0 1129 JMOL=0 1130 do ITYP=1,NTYPES 1131 IBEG=IADDR(ITYP)+1 1132 IEND=IADDR(ITYP+1) 1133 if(IINIT(ITYP).eq.1)then 1134 if(NSPEC(ITYP).gt.1)stop 1135 +' only 1 molecule of such type is allowed' 1136 JMOL=JMOL+1 1137 X(JMOL)=0.d0 1138 Y(JMOL)=0.d0 1139 Z(JMOL)=0.d0 1140 if(IPRINT.ge.8)write(*,*) 1141 +' put molecule ',JMOL,' at zero ref.point' 1142 else 1143 do I=IBEG,IEND 1144 J=J+1 1145 JMOL=JMOL+1 1146 X(JMOL)=XRS(J) 1147 Y(JMOL)=YRS(J) 1148 Z(JMOL)=ZRS(J) 1149 end do 1150 end if 1151 end do 1152 if(TASKID.eq.MAST)write(*,*)JMOL,' molecules are distributed' 1153 if(IPRINT.ge.8)then 1154 WRITE(*,*)' COM of molecules:' 1155 DO ityp=1,ntypes 1156 IBEG=IADDR(ITYP)+1 1157 IEND=IADDR(ITYP+1) 1158 do I=IBEG,IEND 1159 write(*,'(2i5,3f14.5)')ITYP,I,X(I),Y(I),Z(I) 1160 end do 1161 end do 1162 end if 1163 return 1164 end 1165* 1166*=================== SPHOLE ==================================== 1167* 1168* 7. Put one large molecule in a spherical hole and 1169* distribute other molecules around it 1170* 1171 subroutine SPHOLE 1172 include "prcm.h" 1173 dimension XRS(NPART),YRS(NPART),ZRS(NPART) 1174 NOPSM=0 1175 do ITYP=1,NTYPES 1176 if(IINIT(ITYP).ne.1)NOPSM=NOPSM+NSPEC(ITYP) 1177 end do 1178 VOLAV=VOL-4.d0*PI*RDNA**3/3.d0 1179 if(TASKID.eq.MAST) 1180 +write(*,*)' Distribute ',NOPSM,' molecules in ',VOLAV,' A**3' 1181 SHAG=(VOLAV/NOPSM)**0.333333333 1182 IAT=0 1183 NSX=BOXL/SHAG+0.7 1184 NSY=BOYL/SHAG+0.7 1185 NSZ=BOZL/SHAG+0.5 1186 10 NCOUNT=0 1187 do IX=1,NSX 1188 XX=-HBOXL+BOXL*(dfloat(IX)-0.6)/NSX 1189 do IY=1,NSY 1190 YY=-HBOYL+BOYL*(dfloat(IY)-0.45)/NSY 1191 do IZ=1,NSZ 1192 ZZ=-HBOZL+BOZL*(dfloat(IZ)-0.55)/NSZ 1193 RR=sqrt(XX**2+YY**2+ZZ**2) 1194 if(RR.gt.RDNA.and.(.not.LOCT.or. 1195 +abs(XX)+abs(YY)+abs(ZZ).lt.0.75*BOXL))then 1196 NCOUNT=NCOUNT+1 1197 XRS(NCOUNT)=XX 1198 YRS(NCOUNT)=YY 1199 ZRS(NCOUNT)=ZZ 1200 if(NCOUNT.ge.NOPSM)go to 20 1201 end if 1202 end do 1203 end do 1204 end do 1205 if(NCOUNT.le.NOPSM)then 1206 IAT=IAT+1 1207 if(mod(IAT,3).eq.0)then 1208 NSX=NSX+1 1209 if(TASKID.eq.MAST)write(*,*)' increase NSX to ',NSX 1210 else if(mod(IAT,3).eq.1)then 1211 NSY=NSY+1 1212 if(TASKID.eq.MAST)write(*,*)' increase NSY to ',NSY 1213 else 1214 NSZ=NSZ+1 1215 if(TASKID.eq.MAST)write(*,*)' increase NSZ to ',NSZ 1216 end if 1217 go to 10 1218 end if 1219* blandning 1220 20 call MIXTURE(XRS,YRS,ZRS,NOPSM) 1221* put fixed molecules on their places 1222 J=0 1223 JMOL=0 1224 do ITYP=1,NTYPES 1225 IBEG=IADDR(ITYP)+1 1226 IEND=IADDR(ITYP+1) 1227 if(IINIT(ITYP).eq.1)then 1228 if(NSPEC(ITYP).gt.1)stop 1229 +' only 1 molecule of such type is allowed' 1230 JMOL=JMOL+1 1231 X(JMOL)=0.d0 1232 Y(JMOL)=0.d0 1233 Z(JMOL)=0.d0 1234 if(IPRINT.ge.8)write(*,*) 1235 +' put molecule ',JMOL,' at zero ref.point' 1236 else 1237 do I=IBEG,IEND 1238 J=J+1 1239 JMOL=JMOL+1 1240 X(JMOL)=XRS(J) 1241 Y(JMOL)=YRS(J) 1242 Z(JMOL)=ZRS(J) 1243 end do 1244 end if 1245 end do 1246 if(TASKID.eq.MAST)write(*,*)JMOL,' molecules are distributed' 1247 if(IPRINT.ge.8)then 1248 WRITE(*,*)' COM of molecules:' 1249 DO ityp=1,ntypes 1250 IBEG=IADDR(ITYP)+1 1251 IEND=IADDR(ITYP+1) 1252 do I=IBEG,IEND 1253 write(*,'(2i5,3f14.5)')ITYP,I,X(I),Y(I),Z(I) 1254 end do 1255 end do 1256 end if 1257 return 1258 end 1259* 1260*===================== MIXTURE =================================== 1261* 1262* 8. Mix molecules to avoid artificial correlations 1263* 1264 subroutine MIXTURE(X,Y,Z,NOP) 1265 real*8 X(NOP),Y(NOP),Z(NOP) 1266* 1267* making mixture 1268* 1269 do ICT=2,7 1270 if(mod(ICT,2).eq.0)then 1271 J=NOP 1272 do I=1,NOP,ICT 1273 XX=X(I) 1274 YY=Y(I) 1275 ZZ=Z(I) 1276 X(I)=X(J) 1277 Y(I)=Y(J) 1278 Z(I)=Z(J) 1279 X(J)=XX 1280 Y(J)=YY 1281 Z(J)=ZZ 1282 J=J-1 1283 end do 1284 else 1285 J=1 1286 do I=1,NOP,ICT 1287 XX=X(I) 1288 YY=Y(I) 1289 ZZ=Z(I) 1290 X(I)=X(J) 1291 Y(I)=Y(J) 1292 Z(I)=Z(J) 1293 X(J)=XX 1294 Y(J)=YY 1295 Z(J)=ZZ 1296 J=J+1 1297 end do 1298 end if 1299 end do 1300 return 1301 end 1302* 1303*================ PULLINI ======================================= 1304* 1305 subroutine PULLINI 1306C 1307C 9. This subroutine read coordinates of points to which some 1308C atoms will be bound by harmonical forces (see subr PULLBACK) 1309C Normally, it is not used (LRR=.false.) 1310C 1311 include "prcm.h" 1312 character*128 STR,TAKESTR 1313 open(unit=31,file=filref,status='old',err=98) 1314 STR=TAKESTR(31,IE) 1315 read(STR,*,end=99,err=99)NPULL 1316 if(NPULL.gt.NPLM)stop '!!! increase NPLM!!!' 1317 if(IPRINT.ge.6)write(*,*)' Linked atoms: ' 1318 do I=1,NPULL 1319C IS - atom number 1320 STR=TAKESTR(31,IE) 1321 read(STR,*,end=99,err=99)IS,XX,YY,ZZ 1322 INPL(I)=IS 1323 XPL(I)=XX 1324 YPL(I)=YY 1325 ZPL(I)=ZZ 1326 if(IPRINT.ge.6)write(*,'(i5,2x,a,3f14.6)') 1327 + IS,NM(NSITE(IS)),XX,YY,ZZ 1328 end do 1329 close(31) 1330 return 1331 98 if(TASKID.eq.MAST)write(*,*) 1332 +' LRR=.true. and file ',filref,' not found' 1333 call FINAL 1334 99 IE=99 1335 STR=TAKESTR(31,IE) 1336 end 1337* 1338*========== MKPARLIST ======================================= 1339* 1340* Create file with list of pairs for LJ interactions 1341* 1342 subroutine MKPARLIST 1343 include "prcm.h" 1344 character*128 TAKESTR,STR,KEYW*32,CHT1*8,CHT2*8 1345 open(unit=24,file=fpairs,status='old',err=198) 1346 open(unit=25,file='pairs.list',status='unknown') 1347 write(25,'(a)')'# LJ interactions pair list' 1348 IOK=0 1349 10 STR=TAKESTR(24,IE) 1350 if(IE.eq.-1)go to 199 1351 read(STR,*,end=10)KEYW 1352 if(KEYW(1:8).eq.'LJ_PAIRS')then 1353 IOK=1 1354 if(IPRINT.ge.6)write(*,*) 1355 + 'LJ_PAIRS list section found in the force field file' 1356 go to 10 1357 end if 1358 if(KEYW(1:3).eq.'END')then 1359 IOK=2 1360 go to 199 1361 end if 1362 if(IOK.eq.1)then 1363 if(LPRINT.ge.7)write(*,'(a72)')STR(1:72) 1364 read(STR,*,end=197,err=197)CHT1,CHT2,SIG,EPS 1365 do IS=1,NS 1366 do JS=1,NS 1367 if(FFT(IS).eq.CHT1.and.FFT(JS).eq.CHT2)then 1368 write(25,'(2I6,2f12.4)')IS,JS,SIG,EPS 1369 end if 1370 end do 1371 end do 1372 end if 1373 go to 10 1374 197 IE=99 1375 STR=TAKESTR(IO,IE) 1376 198 write(*,*)'!!! Force field file not found: ',FPAIRS 1377 write(*,*) 1378 +'Combination rules will be used for all cross-interactions' 1379 LPAIRS=.false. 1380 199 close(24) 1381 close(25) 1382 open(unit=23,file='pairs.list',status='old') 1383 return 1384 end 1385