1 2!DIR$ FIXED 3* 4* Part of MDynaMix package v.5.2 5* 6* Block TRANAL_BASE - analysis of trajectorisproceedies 7* 8* subroutines: 9* SETTRAJ - prepare list of trajectories and data structures 10* READCONF - reads next configuration 11* READMOL - read information about molecules from .mmol files 12* GETCOM - compute molecular centres of mass 13* PBC - apply periodic boundary conditions 14* 15* This subroutines prepare data structures 16* for subsequent analysis 17* 18 subroutine SETTRAJ 19 include "tranal.h" 20 character*128 FFILTR,PATHDB,FNAME,cdum*8,chr4*4, FMT*20 21 real*4 TX(NTOT),TY(NTOT),TZ(NTOT) 22 save 23 namelist /TRAJ/NFORM,FNAME,PATHDB,NTYPES,NAMOL,NFBEG,NFEND, 24 + ISTEP,BREAKM,LVEL,LFILTR,FFILTR,LXMOL,FXMOL,IPRINT,LBOND,NSPEC, 25 + NSITS,BOXL,BOYL,BOZL,VFAC,LHEX,LOCT 26* requested parameters: 27* 28* NFORM - trajectory format 29* MDYN - binary MDynaMix 30* XMOL - XMOL (xyz) 31* PDBT - PDB trajectory from GROMACS 32* DCDT - DCD binary trajectory from NAMD 33* 34* PATHDB - path to molecular database (.mmol files) 35* NTYPES - number of molecular types 36* NAMOL - names of molecules to analyse 37* NFBEG - number of the first trajectory file 38* NFEND - number of the last trajectory file 39* NSPEC - number of molecules of each type (not necessary in MDYN format) 40* NSITS - number of atoms in a molecule of each type 41* (not necessary if .mmol files are provided)* 42* other things: 43* 44* BOXL,BOYL,BOZL - box sizes. They are used if trajectory file does not 45* contain information on them 46* 47* LHEX - hexagonal cell 48* LOCT - truncated octahedron cell 49* 50* 51* IPRINT - output level 52* BREAKM - which break in trajectory file (in ps ) is allowed 53* 54* LVEL - read velocities from the trajectory files 55* 56* LXMOL - eventually rewrite trajectory to XMOL format 57* FXMOL - base name of XMOL file 58* ISTEP - how often take configurations 59* 60* LFILTR - analysis for particular atoms 61* FFILTR - file name specifying which atoms take for the analysis 62* 63* LBOND - information about bonds is read from .mmol file 64* 65* VFAC - factor to multiply velocities to bring them to A/fs 66* 67* This is default trajectory format 68 NFORM='MDYN' 69 LVEL=.false. 70 LFILTR=.false. 71 LXMOL=.false. 72 LBOND=.false. 73 LHEX=.false. 74 LOCT=.false. 75 ISTEP=1 76 IPRINT=5 77 BREAKM=111. ! ps 78 VFAC=1. 79 LMMOL=.true. 80* Read input data 81 read(*,TRAJ) 82 if(NFORM.eq.'PDBT'.or.NFORM.eq.'DCDT')then 83 if(LVEL)write(*,*)' no velocities in this trajectory format' 84 LVEL=.false. 85 end if 86 if (LHEX)then 87 BOYL=cos30*BOXL 88 BOXY3=0.66666666666666666667d0*BOYL 89 else if(LOCT)then 90 BOYL=BOXL 91 BOZL=BOXL 92 end if 93* defining the length of the file name 94 IL=0 95 do I=1,120 96 ISTR=ichar(FNAME(I:I)) 97 if(ISTR.le.15)then 98 LF=I-1 99 go to 11 100 end if 101 if(FNAME(I:I).ne.' ')IL=I 102 end do 103 LF=IL 104* Prepare for eventual rewriting trajectory files in XMOL format 105 if(LXMOL)then 106 LLNX=0 107 do I=1,128 108 ISTR=ichar(FXMOL(I:I)) 109 if(ISTR.le.15)then 110 LLNX=I-1 111 go to 11 112 end if 113 if(FXMOL(I:I).ne.' ')LLNX=I 114 end do 115 end if 116* Check and reconstruct file list 117 11 continue 118 write(*,*)FNAME(1:LF) 119 NF=NFEND-NFBEG+1 120 J=0 121 do I=1,NF 122 NFL=NFBEG+I-1 123 filnam(I)(1:LF)=FNAME(1:LF) 124 filnam(I)(LF+1:LF+4)='.000' 125 do M=LF+5,128 126 filnam(I)(M:M)=' ' 127 end do 128 if(NFL.lt.10)then 129 write(filnam(I)(LF+4:LF+4),'(i1)')NFL 130 else if(NFL.lt.100)then 131 write(filnam(I)(LF+3:LF+4),'(i2)')NFL 132 else if(NFL.lt.1000)then 133 write(filnam(I)(LF+2:LF+4),'(i3)')NFL 134 else 135 write(filnam(I)(LF+2:LF+5),'(i4)')NFL 136 end if 137 write(*,*)' looking for file ',filnam(I) 138 if(NFORM.eq.'PDBT')then 139 open (unit=12,file=filnam(I),status='old',err=5) 140 if(IPRINT.ge.6)write(*,'(a6,a64,a8)') 141 + ' file ',filnam(I)(1:64),' is open' 142 9 read(12,'(a)',err=10)STR 143* PDB trajectory format 144 if(STR(1:6).ne.'HEADER'.and.STR(1:5).ne.'TITLE')go to 9 145 do K=7,128 146 if(STR(K:K+1).eq.'t=')then 147 read(str(K+2:K+20),*)FULTIM 148 if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps' 149 go to 13 150 end if 151 end do 152 write(*,*)' time is not defined ' 153 13 continue 154 go to 15 155* XMOL trajectory format 156 else if(NFORM.eq.'XMOL')then 157 open (unit=12,file=filnam(I),status='old',err=5) 158 if(IPRINT.ge.6)write(*,'(a6,a64,a8)') 159 + ' file ',filnam(I)(1:64),' is open' 160 read(12,*,err=10)NSTOT 161 read(12,*,err=14)cdum,FULTIM 162 FULTIM=0.001*FULTIM ! in ps 163 if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps' 164 go to 15 165 14 write(*,*)'time is not defined' 166 FULTIM=I 167 go to 15 168* Binary MDynaMix format 169 else if(NFORM.eq.'MDYN')then 170 open (unit=12,file=filnam(I),status='old', 171 + form='unformatted',err=5) 172 if(IPRINT.ge.6)write(*,'(a6,a64,a8)') 173 + ' file ',filnam(I)(1:64),' is open' 174 read(12,err=10,end=10)IA,dt,boxl,boyl,bozl,unitl,ntypes 175 if(IA.eq.1)LVEL=.true. 176 if(IPRINT.ge.5)write(*,*)' num.of types ',NTYPES 177 if(I.gt.1)then 178 if(NTP.ne.NTYPES)then 179 write(*,*) 180 + ' uncompatible NTYPES in ',filnam(I),' and ',filnam(I-1) 181 stop 182 end if 183 else 184 rewind (12) 185 read(12,err=10) IA, dt, boxl, boyl, bozl, unitl, ntypes, 186 x (nspec(ii),nsits(ii),logo,ii=1,ntypes) 187 end if 188 NTP=NTYPES 189 do ITYP=1,NTYPES 190 read(12,err=10) 191 if(IPRINT.ge.8)write(*,*)ITYP,NSPEC(ITYP),NSITS(ITYP) 192 end do 193 read(12)IA,FULTIM 194 FULTIM=FULTIM*1.d12 ! in ps 195 if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps' 196 go to 15 197* NAMD DCD format 198 else if(NFORM.eq.'DCDT')then 199 open (unit=12,file=filnam(I),status='old', 200 + form='unformatted',err=5) 201 if(IPRINT.ge.6)write(*,'(a6,a64,a8)') 202 + ' file ',filnam(I)(1:64),' is open' 203 read(12)CHR4,NCF,NST1,NFREQ,NSTF,I1,I2,I3,I4,I5,DTN 204 write(*,*) CHR4,NCF 205 FULTIM=NST1*DTN*0.04888 ! in ps 206 if(IPRINT.ge.5)write(*,*)' FULTIM ',FULTIM,' ps' 207 go to 15 208 else 209 write(*,*)' Unsupported trajectory format ',NFORM 210 stop 211 end if 212 5 write(*,*)' File not found: ',filnam(I) 213 go to 20 214 10 write(*,*)' File ',filnam(I),': input error' 215 go to 20 216 15 J=J+1 217 TIMIN(J)=FULTIM 218 filnam(J)=filnam(I) 219 20 close(12) 220 end do 221 if(J.eq.0)then 222 write(*,*)' No trajectory files found. ' 223 stop 224 else 225 write(*,*)' input of the file list is completed:',J,' files' 226 end if 227 30 NF=J 228 FULTIM0=TIMIN(1) 229 TTER=TTER*1.d-12 230 if(NSTOT.gt.NTOT)then 231 write(*,*)'total number of atoms in the trajectory ', 232 + NSTOT,' exceds the limit' 233 write(*,*) 234 + ' Increase parameter NTOT in tranal.h file and recompile' 235 end if 236 if(NTYPES.le.0)then 237 write(*,*)'!!! Wrong number of molecule types : ',NTYPES 238 stop 239 end if 240* Adresses and referenses 241 IAN=0 ! number of analysed configurations 242 ISADR(1)=0 243 ISADDR(1)=0 244 IADDR(1)=0 245 NOP=0 246 NSITES=0 247 NSTOT=0 248 FNSTR=0. 249 IATOM=0 250 IMOL=0 251 IST=0 252 NCOR=0 253 JBREAK=1 ! indicate break of trajectory 254 NX=0 255* Read molecular information 256 do ITYP=1,NTYPES 257 LIST(ITYP)=1 258 if(NFORM.eq.'MDYN')NSTS=NSITS(ITYP) 259 CALL READMOL (SUMM,NX,ITYP,PATHDB,NAMOL(ITYP),IER) 260 if(NFORM.eq.'MDYN'.and.NSTS.ne.NSITS(ITYP))then 261 write(*,*)' Molecule of type ',ITYP,' has ',NSTS, 262 &' sites in the trajectory but ',NSITS(ITYP),' in the input file' 263 stop 264 end if 265 NOP = NOP + NSPEC(ITYP) 266 TOTMAS = TOTMAS+SUMM*NSPEC(ITYP) 267 SUMMAS(ITYP) = SUMM 268 NAME(ITYP)=NAMOL(ITYP)(1:6) 269 end do 270 if(NOP.le.0)then 271 write(*,*)'!!! Wrong number of molecules : ',NOP 272 stop 273 end if 274 FACTM=AVSNO*1.d3 275 do ityp=1, ntypes 276 ISADR(ITYP+1) = ISADR(ITYP)+ NSITS(ITYP) ! site addr 277 ISADDR(ITYP+1) = ISADDR(ITYP)+NSPEC(ITYP)*NSITS(ITYP) ! atom adr 278 IADDR (ITYP+1) = IADDR (ITYP)+NSPEC(ITYP) ! mol. addr 279 NSITES = NSITES + NSITS(ITYP) 280 if(NSITS(ITYP).eq.2)FNSTR=FNSTR+2.d0*NSPEC(ITYP) 281 if(NSITS(ITYP).ge.3)FNSTR=FNSTR+3.d0*NSPEC(ITYP) 282 NSTOT = NSTOT + NSPEC(ITYP)*NSITS(ITYP) 283 FNSS3 = 3.D0*DFLOAT(NSITES) 284 write(*,'(a6,i4,3x,a5,i5,4x,a6,i4)') 285 +' type ',ITYP,' Nmol',NSPEC(ITYP), 286 +'Nsites',NSITS(ITYP) 287 do J=ISADR(ITYP)+1,ISADR(ITYP+1) 288 MASSD(J)=MASS(J)/TOTMAS/FACTM 289 ITS(J)=ITYP 290 end do ! Nsite -> type 291 do IS=1,NSPEC(ITYP) 292 IMOL=IMOL+1 293 do I=1,NSITS(ITYP) 294 ISITE=ISADR(ITYP)+I 295 IATOM=IATOM+1 296 NNUM(IATOM)=IMOL ! Natom -> Nmol 297 NSITE(IATOM)=ISITE ! Natom -> Nsite 298 ITYPE(IATOM)=ITYP ! Natom -> type 299 if(IPRINT.ge.9)write(*,'(4I6,2x,2a4,a6)') 300 + IATOM,IMOL,ISITE,ITYP,NM(ISITE) 301 end do 302 ITM(IMOL)=ITYP ! mol -> type 303 end do 304 end do 305* Summary of the system 306 write(*,*)' --------------------------------' 307 write(*,*)' Summary of the system:' 308 if(LMMOL)then 309 write(*,*)' based on data from .mmol files' 310 else 311 write(*,*)' Not all .mmol files were found' 312 write(*,*)' information is generated from the input file' 313* giving numeric site names 314 do IS=1,NSITES 315 MASS(IS)=1. 316 NM(IS)='Art' 317 write(NM(IS),'(i4)')IS 318 end do 319 do ITYP=1,NTYPES 320 SUMMAS(ITYP)=1.*NSITS(ITYP) 321 end do 322 end if 323 write(*,*)' Type Name Atoms Molecules' 324 do ITYP=1,NTYPES 325 write(*,'(i3,3x,a6,2i7)') 326 + ITYP,NAME(ITYP),NSITS(ITYP),NSPEC(ITYP) 327 end do 328 write(*,*)' Total number of atoms: ',NSTOT 329 write(*,*)' --------------------------------' 330* Renumeration array 331 if(LFILTR)then 332 do I=1,NUMTR 333 NUMR(I)=0 334 end do 335 write(*,*)' renumerating atoms from file ',FFILTR 336 open(unit=14,file=FFILTR,status='old') 337 do I=1,NSTOT 338 read(14,*,err=921,end=921)II,NUMO(II) 339 NUMR(NUMO(II))=II 340 if(IPRINT.ge.7)write(*,*)I,II,NUMO(II) 341 end do 342 close(14) 343 else 344 NUMTR=NSTOT 345 do I=1,NSTOT 346 NUMO(I)=I 347 NUMR(I)=I 348 end do 349 end if 350 FULTIM1=TIMIN(1) 351 write(*,*)' FULTIM1=',FULTIM1 352 return 353 921 write(*,*)' error reading renumeration file ',FFILTR 354 stop 355 end 356* 357*============= READCONF =========================== 358* 359* read configuration 360* IEND signals the end of trajectories 361 subroutine READCONF(IEND) 362 include "tranal.h" 363 real*4 TX(NTOT),TY(NTOT),TZ(NTOT) 364 character cdum*32,chr4*4 365 data IINIT/0/ 366 save 367 if(IINIT.eq.0)then 368 IINIT=1 369 IFL=1 ! number of the first file to read 370 IIN=0 ! Signals that we read headers 371 end if 372 IEND=0 ! signals end of the analysis 373 ITRI=0 ! signals about double point 374* Begin reading of the next trajectory file 375* Reading HEADERs 376* Case of PDB trajactory 377 35 if(NFORM.eq.'PDBT')then 378 if(IIN.eq.0)then 379 open (unit=12,file=filnam(IFL),status='old', 380 + form='formatted') 381 if(IPRINT.ge.7)write(*,*)' file ',IFL,' open' 382 end if 383 31 read(12,'(a)',end=45) STR 384 if(STR(1:6).ne.'HEADER'.and.STR(1:5).ne.'TITLE')go to 31 385 do K=7,128 386 if(STR(K:K+1).eq.'t=')then 387 read(str(K+2:K+20),*)FULTIM 388 go to 33 389 end if 390 end do 391 write(*,*)' time is not defined ' 392 33 IFLR=0 393* Line 2 394 32 read( 12,'(a)',end=45 ) STR 395 if(STR(1:4).eq.'ATOM')then 396 if(IPRINT.ge.6)write(*,*)' box sizes undetermined' 397 IFLR=1 398 go to 34 399 end if 400 if(STR(1:4).ne.'CRYS')go to 32 401 read(str,*)cdum,BOXL,BOYL,BOZL 402 if(IPRINT.ge.7)write(*,'(a,3f12.3)')' Box: ',BOXL,BOYL,BOZL 403 else if(NFORM.eq.'XMOL')then 404* XMOL trajectory 405 if(IIN.eq.0)then 406 open (unit=12,file=filnam(IFL),status='old', 407 + form='formatted') 408 if(IPRINT.ge.7)write(*,*)' file ',IFL,' open' 409 end if 410 read(12,*,err=40,end=45) NUMTR 411 read( 12,'(a)',err=40,end=40 )STR 412 read(STR,*,end=40,err=36)cdum,FULTIM 413 do J=1,100 414 if(STR(J:J+2).eq.'BOX')then 415 read(STR(J+4:128),*,end=40,err=40)BOXL,BOYL,BOZL 416 if(IPRINT.ge.7)write(*,'(a,3f12.3)') 417 & ' Box: ',BOXL,BOYL,BOZL 418 go to 44 419 end if 420 end do 421 if(IPRINT.ge.6)write(*,*)' box sizes undetermined' 422 44 do J=20,110 423 if(STR(J:J+4).eq.'EEnum')then 424 read(STR(J+5:128),*,end=40,err=40)MEE 425 if(IPRINT.ge.7)write(*,'(a,3f12.3)') 426 & ' EE: ',MEE 427 go to 37 428 end if 429 end do 430 go to 37 431 36 if(IPRINT.ge.6)write(*,*)' time is undetermined' 432 FULTIM=IAN 433 37 FULTIM=0.001*FULTIM ! in ps 434 else if(NFORM.eq.'DCDT')then 435* NAMD DCD trajectory 436 if(IIN.eq.0)then 437 open (unit=12,file=filnam(IFL),status='old', 438 + form='unformatted') 439 read(12,err=40,end=45) 440 + CHR4,NCF,NST1,NFREQ,NSTF,I1,I2,I3,I4,I5,DTN 441 FULTIM=NST1*DTN*0.04888 ! in ps 442 read(12,err=40,end=45) 443 read(12,err=40,end=45)NUMTR 444 else 445 FULTIM=FULTIM+NFREQ*DTN*0.04888 446 end if 447 read(12,err=40,end=45)BOXL,Y1,BOYL,Y2,Y3,BOZL 448 if(IPRINT.ge.7)write(*,'(a,3f12.3)')' Box: ',BOXL,BOYL,BOZL 449 if(LXMOL)then 450 filxm(1:LLNX)=FXMOL(1:LLNX) 451 filxm(LLNX+1:LLNX+4)=filnam(IFL)(LF+1:LF+4) 452 open(unit=19,file=filxm,status='unknown') 453 end if 454 else 455* Binary MDynamix trajectory 456 if(IIN.eq.0)then 457 open (unit=12,file=filnam(IFL),status='old', 458 + form='unformatted',err=45) 459 read(12) IVEL, dt, boxl, boyl, bozl, unitl, ntypes, 460 x (nspec(i),nsits(i),logo,i=1,ntypes) 461 if(LXMOL)then 462 filxm(1:LLNX)=FXMOL(1:LLNX) 463 filxm(LLNX+1:LLNX+4)=filnam(IFL)(LF+1:LF+4) 464 open(unit=19,file=filxm,status='unknown') 465 end if 466 if(IPRINT.ge.7)write(*,*) 467 & 'read first line.',NTYPES,' types' 468 do ityp=1, ntypes 469 NSB = ISADR (ITYP)+1 470 NSE = ISADR (ITYP+1) 471* Lines from 2 to NTYPES+1 472 read( 12 ) (mass(j), j=NSB,NSE),LIST(ITYP) 473* LIST(ITYP)=1 474 if(IPRINT.ge.8)write(*,*)ityp,NSB,NSE,LIST(ITYP) 475 end do 476 if(IPRINT.ge.7)write(*,*)' read masses ' 477 end if 478* First line of a configureation 479 read(12,err=40,end=45)MEE,fultim,cumtim, 480 & TEMP,PRES,EP,BOXL,BOYL,BOZL,(LIST(LL),LL=1,NTYPES) 481 FULTIM=FULTIM*1.d12 ! in ps 482 if(IPRINT.ge.7)write(*,'(2(a,3f10.2),a,I3)') 483 & 'Box: ',BOXL,BOYL,BOZL,' term:',TEMP,PRES,EP,' EE ',MEE 484 end if ! if(NFORM... 48534 if(IIN.eq.0)then 486 write(*,'(2a/a,f14.3,a8,6I3)') 487 + ' begin analys of file ',filnam(IFL), 488 + ' init. T = ',fultim 489 IIN=1 490 end if 491 IBREAK=0 492 DIFT=FULTIM-FULTIM1 493 HBOXL=0.5*BOXL 494 HBOYL=0.5*BOYL 495 HBOZL=0.5*BOZL 496 VOL=BOXL*BOYL*BOZL 497 if(DIFT.gt.BREAKM)then 498 write(*,*) 499 + 'break in trajectory file. from ',FULTIM1,' to ', 500 + FULTIM,'ps' 501 IBREAK=1 502 end if 503 if(IFL.lt.NF.and.fultim.ge.TIMIN(IFL+1))go to 45 504 FULTIM1=FULTIM 505 if(IPRINT.ge.7)write(*,*)' reading conf at time ',FULTIM 506 if(NFORM.eq.'MDYN')then 507c - for each type 508 do ityp=1,ntypes 509 NSB = ISADDR(ITYP)+1 510 NSE = ISADDR(ITYP+1) 511* Other NTYPES lines of each point 512 if(LIST(ITYP).ne.0)then 513 if(IVEL.le.0)then 514 read(12,err=40,end=45) (TX(j),j=NSB,NSE), 515 * (TY(j),j=NSB,NSE), 516 * (TZ(j),j=NSB,NSE) 517 do J=NSB,NSE 518 SX(J)=TX(J) 519 SY(J)=TY(J) 520 SZ(J)=TZ(J) 521 end do 522 else ! IVEL=1 523 read(12,err=40,end=45) (TX(j),j=NSB,NSE), 524 * (TY(j),j=NSB,NSE), 525 * (TZ(j),j=NSB,NSE) 526 do J=NSB,NSE 527 SX(J)=TX(J) 528 SY(J)=TY(J) 529 SZ(J)=TZ(J) 530 end do 531 read(12,err=40,end=45) (TX(j),j=NSB,NSE), 532 * (TY(j),j=NSB,NSE), 533 * (TZ(j),j=NSB,NSE) 534 do J=NSB,NSE 535 VX(J)=TX(J) 536 VY(J)=TY(J) 537 VZ(J)=TZ(J) 538 end do 539 end if ! IVEL 540 end if ! LIST 541 end do ! ITYP 542 else if(NFORM.eq.'DCDT')then 543 read(12,err=40,end=45)(TX(I),I=1,NUMTR) 544 read(12,err=40,end=45)(TY(I),I=1,NUMTR) 545 read(12,err=40,end=45)(TZ(I),I=1,NUMTR) 546 do i=1,NUMTR 547 OX(I)=SX(I) 548 OY(I)=SY(I) 549 OZ(I)=SZ(I) 550 SX(I)=TX(I) 551 SY(I)=TY(I) 552 SZ(I)=TZ(I) 553 end do 554 else ! ASCII trajectories 555 do i=1,NUMTR 556 OX(I)=SX(I) 557 OY(I)=SY(I) 558 OZ(I)=SZ(I) 559 IN=NUMR(I) 560 if(LVEL.and.NFORM.eq.'XMOL')then 561 read(12,*,err=40,end=40)CHR4,XX,YY,ZZ,VVX,VVY,VVZ 562 if(.not.LMMOL)NM(NSITE(IN))=CHR4 563 if(IN.ne.0)then 564 SX(IN)=XX 565 SY(IN)=YY 566 SZ(IN)=ZZ 567 VX(IN)=VVX 568 VY(IN)=VVY 569 VZ(IN)=VVZ 570 if(IPRINT.ge.8)write(*,'(2i5,3f10.4)')IN,I,XX,YY,ZZ 571 end if 572 else 573* GROMACS PDB / XMOL trajectory 574 if(NFORM.eq.'PDBT')then 575 53 read(12,'(a)',err=40)STR 576 if(STR(1:6).eq.'ENDMDL'.or.STR(1:6).eq.'HEADER')then 577 write(*,*)' Number of atoms in the trajectory ',IN-1, 578 & ' < than in the input file ',NSTOT 579 stop 580 end if 581 if(STR(1:4).ne.'ATOM')go to 53 582 read(str(31:128),*)XX,YY,ZZ 583 if(.not.LMMOL)NM(NSITE(IN))(1:3)=STR(14:16) 584 else ! XMOL case 585 read(12,*,err=40,end=40)CHR4,XX,YY,ZZ 586 if(.not.LMMOL)NM(NSITE(IN))=CHR4 587 end if 588 if(IN.ne.0)then 589 SX(IN)=XX 590 SY(IN)=YY 591 SZ(IN)=ZZ 592 if(IPRINT.ge.8)write(*,'(2i5,3f10.4)')IN,I,XX,YY,ZZ 593 end if 594 end if 595 end do 596 end if 597 JBREAK=0 598 if(ITRI.ne.0.and.DIFT.le.0.00001)then 599 write(*,*)' double point in the trajectory file at ',FULTIM 600 go to 35 601 end if 602 ITRI=1 603 DT=DIFT 604 if(DIFT.le.0.0001)DIFT=0.0001 605 EKIN=0. 606 if(LOCT) VOL=0.5*VOL 607 do I=1,NSTOT 608 DX=SX(I)-OX(I) 609 DY=SY(I)-OY(I) 610 DZ=SZ(I)-OZ(I) 611 if(DX.gt.HBOXL)DX=DX-BOXL 612 if(DX.lt.-HBOXL)DX=DX+BOXL 613 if(DY.gt.HBOYL)DY=DY-BOZL 614 if(DY.lt.-HBOYL)DY=DY+BOZL 615 if(DZ.gt.HBOZL)DZ=DZ-BOZL 616 if(DZ.lt.-HBOZL)DZ=DZ+BOZL 617 IS=NSITE(I) 618 if(.not.LVEL)then 619 VX(I)=DX*VFAC*1.d-3/DIFT ! in A/fs 620 VY(I)=DY*VFAC*1.d-3/DIFT 621 VZ(I)=DZ*VFAC*1.d-3/DIFT 622 end if 623 EKIN=EKIN+MASS(IS)*(VX(I)**2+VY(I)**2+VZ(I)**2) ! 2 Ekin 624 OX(I)=SX(I) 625 OY(I)=SY(I) 626 OZ(I)=SZ(I) 627 end do 628 if(LVEL)TEMP=EKIN/(3.d-7*AVSNO*NSTOT*BOLTZ) 629 IEND=0 630 if(mod(IAN,ISTEP).eq.0.and.LXMOL)then 631 NSAT=0 632 do ITYP=1,NTYPES 633 if(LIST(ITYP).ne.0)NSAT=NSAT+NSITS(ITYP)*NSPEC(ITYP) 634 end do 635 write(19,*)NSAT 636 write(19,'(a,f15.2,a,3f11.5)') 637 +' after ',FULTIM*1.d3,' fs, BOX:',BOXL,BOYL,BOZL 638 do I=1,NSTOT 639 ITYP=ITYPE(I) 640 IS=NSITE(I) 641 if(LIST(ITYP).ne.0)write(19,'(a2,3(1x,f9.4))') 642 + NM(IS)(1:2),SX(I),SY(I),SZ(I) 643 end do 644 end if 645 IAN=IAN+1 646 if(BOXL.le.0..or.BOYL.le.0..or.BOZL.le.0.)write(*,*) 647 +'!!! problem with BOX size:',BOXL,BOYL,BOZL 648 return 649* signal reading error 650 40 write(*,*)' input error in file ',filnam(IFL),' T=',FULTIM 651 45 write(*,'(2a/a,f10.3,a,i10)') 652 +' finish analys of file ',filnam(IFL), 653 + ' final T = ',fultim,' N conf. ',IAN 654 close(12) 655 IFL=IFL+1 656 IIN=0 657 if(IFL.gt.NF)then 658 IEND=1 659 return 660 end if 661 go to 35 662 end 663* 664*================ READMOL ============================================= 665* 666* Define molecule from DATABASE 667* 668 subroutine READMOL(SUMM,NX,NTYP,PATHDB,NAMN,IER) 669 include "tranal.h" 670* 671 DIMENSION CM(3) 672 character PATHDB*128,namn*128,FIL*256,FMT*20 673 FIL=' ' 674 NX0=NX 675 IER = 0 676 NRBB = 0 677 IE=0 678 IL=0 679 do I=1,128 680 ISTR=ichar(PATHDB(I:I)) 681 if(ISTR.le.15)then 682 LD=I-1 683 go to 11 684 end if 685 if(PATHDB(I:I).ne.' ')IL=I 686 end do 687 LD=IL 688 11 continue 689 FIL(1:LD)=PATHDB(1:LD) 690 FIL(LD+1:LD+1)='/' 691 IL=0 692 do I=1,32 693 ISTR=ichar(NAMN(I:I)) 694 if(ISTR.le.15)then 695 LN=I-1 696 go to 11 697 end if 698 if(NAMN(I:I).ne.' ')IL=I 699 end do 700 LN=IL 701 12 continue 702 FIL(LD+2:LD+LN+1)=NAMN(1:LN) 703 FIL(LD+LN+2:LD+LN+6)='.mmol' 704 open(unit=12,file=fil,status='old',err=999) 705 if(IPRINT.ge.5)then 706 PRINT "(80('-'))" 707 write(*,*) 708 write(*,'(a,i4,a,a)') 709 + '*** MOLECULE OF TYPE No.',NTYP,' ',NAMN(1:LN) 710 end if 711 STR=TAKESTR(12,IE) 712 read(STR,*)NSTS 713 if(IPRINT.ge.6)PRINT * ,'*** NR OF SITES: ',NSTS 714 if(NX+NSTS.gt.NS) 715 & stop '!!! Increase NS (number of sites) in tranal.h' 716 NSITS(NTYP)=NSTS 717 NSBEG=NX+1 718 do I=1,NSTS 719 STR=TAKESTR(12,IE) 720 IX=NX+I 721! NM(IX)=STR(1:len(NM(IX))) 722! read(STR((len(NM(IX))+1):128),*,err=903)(R(IX,J),J=1,3), 723 read(STR(1:128),*,err=903)NM(IX),(R(IX,J),J=1,3), 724 + MASS(IX),CHARGE(IX) 725 if(IPRINT.ge.6) then 726 write(FMT,*) len(NM(IX)) 727 write(*,'(a5,a'//ADJUSTL(FMT)//',a4,3f7.3,a3,f8.4,a3,f6.3, 728 + a5,f8.4,a4,f9.4)') 729 + 'atom ',NM(IX),' R=',(R(IX,J),J=1,3), 730 + ' M=',MASS(IX),' Q=',CHARGE(IX) 731 732! + write(*,'(a5,a<len(NM(IX))>,a4,3f7.3,a3,f8.4,a3,f6.3, 733! + a5,f8.4,a4,f9.4)') 734! + 'atom ',NM(IX),' R=',(R(IX,J),J=1,3), 735! + ' M=',MASS(IX),' Q=',CHARGE(IX) 736 endif 737! read(STR(5:128),*,err=903)NM(IX),(R(IX,J),J=1,3), 738! + MASS(IX),CHARGE(IX) 739! if(IPRINT.ge.6) 740! + write(*,'(a5,a,a4,3f7.3,a3,f8.4,a3,f6.3,a5,f8.4,a4,f9.4)') 741! + 'atom ',NM(IX),' R=',(R(IX,J),J=1,3), 742! + ' M=',MASS(IX),' Q=',CHARGE(IX) 743 end do 744 NX=NX+NSTS 745 NSEND=NX 746* Reading eventual bond list 747 if(LBOND)then 748 STR=TAKESTR(12,IE) 749 read(STR,*,err=904,end=904)NSR 750 do I=1,NSR 751 STR=TAKESTR(12,IE) 752 if(IPRINT.ge.6)write(*,'(a128)')STR 753 end do 754 STR=TAKESTR(12,IE) 755C Number of bonds 756 read(STR,*,err=902,end=902)NBD 757 NRB(NTYP)=NBD 758 do K=NRBB+1,NRBB+NBD 759 STR=TAKESTR(12,IE) 760 read(STR,*,err=902,end=902)IDUM,IB(K),JB(K) 761 end do 762 NRBB = NRBB+NBD 763 IADB(NTYP)=NRBB-NRB(NTYP) 764 end if 765 close(12) 766* COM coordinates 767 SUMM = 0.D0 768 DO IS = NSBEG,NSEND 769 SUMM = SUMM+MASS(IS) 770 END DO! OF IS 771 SUMMAS(NTYP)=SUMM 772* 773 CM(1) = 0.D0 774 CM(2) = 0.D0 775 CM(3) = 0.D0 776 DO IS = NSBEG,NSEND 777 CM(1) = CM(1)+R(IS,1)*MASS(IS) 778 CM(2) = CM(2)+R(IS,2)*MASS(IS) 779 CM(3) = CM(3)+R(IS,3)*MASS(IS) 780 END DO! OF IS 781* 782 CM(1) = CM(1)/SUMM 783 CM(2) = CM(2)/SUMM 784 CM(3) = CM(3)/SUMM 785 if(IPRINT.ge.10)then 786 PRINT *,'*** Molecular GEOMETRY (C.O.M. IN ORIGIN): ' 787 DO IS = NSBEG,NSEND 788 PRINT "('*** ',I4,3x,A4,3X,3(1X,F7.3))", 789 + IS,NM(IS),R(IS,1),R(IS,2),R(IS,3) 790 END DO! OF IS 791 end if 792 return 793 903 write(*,*)' Cannot read .mmol file for molecule ',NAMN 794 write(*,*)' Error in the list of atoms, atom ',I 795 stop 796 902 write(*,*)' Cannot read .mmol file for molecule ',NAMN 797 write(*,*)' Error in the list of bonds, line ',K 798 stop 799 904 write(*,*)' Cannot read .mmol file for molecule ',NAMN 800 write(*,*)' Error in the reference of .mol file, line ' 801 stop 802 999 write(*,*)'!!! Molecule ',NAMN,' not found in the data base' 803 write(*,*)' File not found: ',FIL 804 write(*,*) 805 +' Will try to retrieve necessary information from other data' 806 LMMOL=.false. 807 return 808 end 809* 810*============== TAKESTR ========================================= 811* 812 character*128 function TAKESTR(KAN,IE) 813 character*128 AUX, fn*32 814 integer LCOUNT(256) 815 data LCOUNT /256*0/ 816 save AUX,LCOUNT 817 if(IE.eq.99)go to 10 818 1 LCOUNT(KAN)=LCOUNT(KAN)+1 819 read(KAN,'(a128)',err=10,end=20)AUX 820 if((AUX(1:1).eq.'#').or.(AUX(1:1).eq.'!'))go to 1 821 TAKESTR=AUX 822 IE=0 823 return 824 10 write(*,*)'!!! error in input file ' 825 if(KAN.eq.5)then 826 write(*,*)' Standard input ' 827 else 828 inquire(unit=KAN,name=fn) 829 write(*,*)' File ',fn 830 end if 831 write(*,*)' in line ',LCOUNT(KAN) 832 write(*,*)AUX 833 stop 834 20 if(IE.ge.0)then 835 write(*,*)'!!! end of file reached' 836 if(KAN.eq.5)then 837 write(*,*)' Standard input ' 838 else 839 inquire(unit=KAN,name=fn) 840 write(*,*)' File ',fn 841 end if 842 stop 843 end if 844 TAKESTR=' ' 845 IE=-1 846 return 847 end 848* 849*===================== UTILITIES ================================= 850* 851* 852*=============== GETCOM ============================================ 853* 854* Compute molecular centres of mass 855* -------------------------------------------------------------------- 856C input: SX,SY,SZ coordinates 857C output: X,Y,Z - molecular COMc 858C WX,WY,WZ - atom coordinates relative to molecular COMs 859C PX,PY,PZ - molecular momenta 860C QX,QY,QZ - molecular angular momenta 861C 862 SUBROUTINE GETCOM 863 include "tranal.h" 864* 865 N = 0 866 I = 0 867 DO ITYP = 1,NTYPES ! over types 868 NSBEG = ISADR (ITYP)+1 869 NSEND = ISADR (ITYP +1) 870 SUMM = SUMMAS (ITYP) 871 DO J = 1,NSPEC(ITYP) ! over molecules 872 I = I+1 873 X (I) = 0.D0 874 Y (I) = 0.D0 875 Z (I) = 0.D0 876 PX (I) = 0.D0 877 PY (I) = 0.D0 878 PZ (I) = 0.D0 879C Calculate C.O.M. vectors relative to the first atom 880 N0 = N+1 881 DO IS = NSBEG,NSEND ! over sites 882 N = N+1 883 DX =SX(N)-SX(N0) 884 DY =SY(N)-SY(N0) 885 DZ =SZ(N)-SZ(N0) 886 call PBC(DX,DY,DZ) 887 X (I) = X(I)+MASS(IS)*DX 888 Y (I) = Y(I)+MASS(IS)*DY 889 Z (I) = Z(I)+MASS(IS)*DZ 890 PX (I) = PX(I)+VX(N)*MASS(IS) 891 PY (I) = PY(I)+VY(N)*MASS(IS) 892 PZ (I) = PZ(I)+VZ(N)*MASS(IS) 893CD if(NNUM(N).ne.I) 894CD +write(*,*)' wrong atom/mol in GETCOM -1: ',N,NNUM(N),I,IS,TASKID 895 END DO! OF IS 896 X (I) = X(I)/SUMM+SX(N0) 897 Y (I) = Y(I)/SUMM+SY(N0) 898 Z (I) = Z(I)/SUMM+SZ(N0) 899 call PBC(X(I),Y(I),Z(I)) 900 END DO! OF I 901 END DO! OF ITYP 902* 903* 7.2 Calculate local atom coordinates relative to molecular COM: 904* -------------------------------------------------------------- 905 N = 0 906 I = 0 907 DO ITYP = 1,NTYPES 908 NSBEG = ISADR (ITYP)+1 909 NSEND = ISADR (ITYP +1) 910 DO J = 1,NSPEC(ITYP) 911 I = I+1 912 QX(I) = 0.D0 913 QY(I) = 0.D0 914 QZ(I) = 0.D0 915 DO IS = NSBEG,NSEND 916 N = N+1 917Calculate short site vectors 918 WX(N) = SX(N)-X(I) 919 WY(N) = SY(N)-Y(I) 920 WZ(N) = SZ(N)-Z(I) 921 call PBC(WX(N),WY(N),WZ(N)) 922 QX(I) = QX(I)+(VY(N)*WZ(N)-WY(N)*VZ(N))*MASS(IS) 923 QY(I) = QY(I)+(VZ(N)*WX(N)-WZ(N)*VX(N))*MASS(IS) 924 QZ(I) = QZ(I)+(VX(N)*WY(N)-WX(N)*VY(N))*MASS(IS) 925CD if(NNUM(N).ne.I) 926CD +write(*,*)' wrong atom/mol in GETCOM -2: ',N,NNUM(N),I,IS,TASKID 927 END DO! OF IS 928 END DO! OF I 929 END DO! OF ITYP 930* reset molecular to COM 931* 932 do I=1,NSTOT 933 IMOL=NNUM(I) 934 SX(I)=X(IMOL)+WX(I) 935 SY(I)=Y(IMOL)+WY(I) 936 SZ(I)=Z(IMOL)+WZ(I) 937 end do 938 RETURN 939 END 940* 941*======================= PBC ======================================== 942* 943 subroutine PBC(XX,YY,ZZ) 944 include "tranal.h" 945 if(XX.gt. HBOXL)XX=XX-BOXL 946 if(XX.lt.-HBOXL)XX=XX+BOXL 947 if(LHEX)then 948C hexagonal periodic cell - minimum image up to 1.5*BOXL from the center 949 XY=BOXYC*XX 950 if(YY.gt.BOXY3-XY.and.(XX.gt.0..or.YY.gt.2.*BOXY3-XY))then 951 YY=YY-BOYL 952 XX=XX-HBOXL 953 XY=BOXYC*XX 954 end if 955 if(YY.lt.-BOXY3-XY.and.(XX.lt.0..or.YY.lt.-2.*BOXY3-XY))then 956 YY=YY+BOYL 957 XX=XX+HBOXL 958 XY=BOXYC*XX 959 end if 960 if(YY.gt.BOXY3+XY)then 961 YY=YY-BOYL 962 XX=XX+HBOXL 963 XY=BOXYC*XX 964 end if 965 if(YY.lt.-BOXY3+XY)then 966 YY=YY+BOYL 967 XX=XX-HBOXL 968 end if 969 else 970C rectangular cell 971 if(YY.gt. HBOYL)YY=YY-BOYL 972 if(YY.lt.-HBOYL)YY=YY+BOYL 973 end if 974 if(ZZ.gt. HBOZL)ZZ=ZZ-BOZL 975 if(ZZ.lt.-HBOZL)ZZ=ZZ+BOZL 976 if(.not.LOCT)return 977C truncated octahedron abs(x)+abs(y)+abs(z) < 0.75*BOXL 978 CORSS=HBOXL*int((4./3.)*(abs(XX)+abs(YY)+abs(ZZ))/BOXL) 979 XX=XX-sign(CORSS,XX) 980 YY=YY-sign(CORSS,YY) 981 ZZ=ZZ-sign(CORSS,ZZ) 982 return 983 end 984* 985*================ LENG =============================================== 986* 987 function LENG(STR,LS) 988 character*1 STR(LS) 989 IL=0 990 do I=1,LS 991 ISTR=ichar(STR(I)) 992 if(ISTR.le.15)then 993 LENG=I-1 994 return 995 end if 996 if(STR(I).ne.' ')IL=I 997 end do 998 LENG=IL 999 return 1000 end 1001* 1002* 1003*============= VECT ========================================= 1004* 1005 subroutine VECT(R1,R2,R3) 1006 real*8 R1(3),R2(3),R3(3) 1007 R3(1)= R1(2)*R2(3)-R1(3)*R2(2) 1008 R3(2)=-R1(1)*R2(3)+R1(3)*R2(1) 1009 R3(3)= R1(1)*R2(2)-R1(2)*R2(1) 1010 return 1011 end 1012* 1013*================ NORM ======================================= 1014* 1015 subroutine NORM(A,B,R,N) 1016 real*8 R,A(*),B(*) 1017 R=0.d0 1018 do I=1,N 1019 R=R+A(I)**2 1020 end do 1021 R=sqrt(R) 1022 if(R.eq.0.d0)then 1023 do I=1,N 1024 B(I)=0.d0 1025 end do 1026 else 1027 do I=1,N 1028 B(I)=A(I)/R 1029 end do 1030 end if 1031 return 1032 end 1033*.... 1034*======== XINTGR ======================================================= 1035*.... 1036 FUNCTION XINTGR(H,Y,NDIM) 1037 IMPLICIT REAL*8 (A-H,O-Z) 1038 DIMENSION Y(0:NDIM) 1039 HT=H/3.0 1040 IF(NDIM.LT.5)RETURN 1041 SUM1=4.00*Y(1) 1042 SUM1=HT*(Y(0)+SUM1+Y(2)) 1043 AUX1=4.00*Y(3) 1044 AUX1=SUM1+HT*(Y(2)+AUX1+Y(4)) 1045 AUX2=HT*(Y(0)+3.8750*(Y(1)+Y(4))+2.6250*(Y(2)+Y(3))+Y(5)) 1046 SUM2=4.00*Y(4) 1047 SUM2=AUX2-HT*(Y(3)+SUM2+Y(5)) 1048 AUX=4.00*Y(2) 1049 IF(NDIM-6)5,5,2 1050 2 DO 4 I=6,NDIM,2 1051 SUM1=AUX1 1052 SUM2=AUX2 1053 AUX1=4.00*Y(I-1) 1054 AUX1=SUM1+HT*(Y(I-2)+AUX1+Y(I)) 1055 IF(I-NDIM)3,6,6 1056 3 AUX2=4.00*Y(I) 1057 AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1)) 1058 4 CONTINUE 1059 5 CONTINUE 1060 XINTGR=AUX2 1061 RETURN 1062 6 CONTINUE 1063 XINTGR=AUX1 1064 RETURN 1065 END 1066* 1067*===================== INV3B3 ========================================== 1068* 1069* Invert matrix 3x3 1070* 1071 SUBROUTINE INV3B3(A,AI,DET,IER) 1072 IMPLICIT real*8 (A-H,O-Z) 1073 DIMENSION A(9),AI(9) 1074* 1075 AI(1) = A(5)*A(9)-A(6)*A(8) 1076 AI(2) = A(3)*A(8)-A(2)*A(9) 1077 AI(3) = A(2)*A(6)-A(3)*A(5) 1078 AI(4) = A(6)*A(7)-A(4)*A(9) 1079 AI(5) = A(1)*A(9)-A(3)*A(7) 1080 AI(6) = A(3)*A(4)-A(1)*A(6) 1081 AI(7) = A(4)*A(8)-A(5)*A(7) 1082 AI(8) = A(2)*A(7)-A(1)*A(8) 1083 AI(9) = A(1)*A(5)-A(2)*A(4) 1084* 1085 DET = A(1)*AI(1)+A(4)*AI(2)+A(7)*AI(3) 1086 SPUR=A(1)+A(5)+A(9) 1087 IER=0 1088 if(dabs(DET).lt.1.d-6*dabs(SPUR**3))IER=1 1089 IF(DABS(DET).gt.0.D0) then 1090 R=1.D0/DET 1091 else 1092 R=0.d0 1093 IER=2 1094 end if 1095* 1096 AI(1) = R*AI(1) 1097 AI(2) = R*AI(2) 1098 AI(3) = R*AI(3) 1099 AI(4) = R*AI(4) 1100 AI(5) = R*AI(5) 1101 AI(6) = R*AI(6) 1102 AI(7) = R*AI(7) 1103 AI(8) = R*AI(8) 1104 AI(9) = R*AI(9) 1105* 1106 RETURN 1107 END 1108