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