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