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