/*************************************************************************/
/*                                                                       */
/*  Copyright (c) 1994 Stanford University                               */
/*                                                                       */
/*  All rights reserved.                                                 */
/*                                                                       */
/*  Permission is given to use, copy, and modify this software for any   */
/*  non-commercial purpose as long as this copyright notice is not       */
/*  removed.  All other uses, including redistribution in whole or in    */
/*  part, are forbidden without prior written permission.                */
/*                                                                       */
/*  This software is provided with absolutely no warranty and no         */
/*  support.                                                             */
/*                                                                       */
/*************************************************************************/

EXTERN_ENV
#include "math.h"
#include "mdvar.h"
#include "water.h"
#include "wwpot.h"
#include "cnst.h"
#include "parameters.h"
#include "mddata.h"
#include "split.h"
#include "global.h"

extern double ****PFORCES;

/* in this version of interf, a private force array is maintained */
/* for every process.  A process computes interactions into its   */
/* private force array, and later updates the shared destination  */
/* array with locks, but only updates those locations that it     */
/* computed something for                                         */

INTERF(DEST,VIR,ProcID)
  int DEST;
  double *VIR;
  unsigned ProcID;
{
    /* This routine gets called both from main() and from mdmain().
       When called from main(), it is used to estimate the initial
       accelerations by computing intermolecular forces.  When called
       from mdmain(), it is used to compute intermolecular forces.
       The parameter DEST specifies whether results go into the 
       accelerations or the forces. Uses routine UPDATE_FORCES in this
       file, and routine CSHIFT in file cshift.U */
    /*
      .....this routine calculates inter-molecular interaction forces
      the distances are arranged in the order  M-M, M-H1, M-H3, H1-M,
      H3-M, H1-H3, H1-H1, H3-H1, H3-H3, O-O, O-H1, O-H3, H1-O, H3-O, 
      where the M are "centers" of the molecules.
      */
    
    int mol, comp, dir, icomp;
    int comp_last, half_mol;
    int    KC, K;
    double YL[15], XL[15], ZL[15], RS[15], FF[15], RL[15]; /* per-
                                                              interaction arrays that hold some computed distances */
    double  FTEMP;
    double LVIR = 0.0;
    double *temp_p;
    
    { /* initialize PFORCES array */
        
        int ct1,ct2,ct3;
        
        for (ct1 = 0; ct1<NMOL; ct1++)
            for (ct2 = 0; ct2<NDIR; ct2++)
                for (ct3 = 0; ct3<NATOM; ct3++) 
                    PFORCES[ProcID][ct1][ct2][ct3] = 0; 
        
    } 
    
    half_mol = NMOL/2;
    for (mol = StartMol[ProcID]; mol < StartMol[ProcID+1]; mol++) {
        comp_last = mol + half_mol;
        if (NMOL%2 == 0) {
            if ((half_mol <= mol) && (mol%2 == 0)) { 
                comp_last--;
            }
            if ((mol < half_mol) && (comp_last%2 == 1)) {
                comp_last--;
            }
        }
        for (icomp = mol+1; icomp <= comp_last; icomp++) {
            comp = icomp;
            if (comp > NMOL1) comp = comp%NMOL;
            
            /*  compute some intermolecular distances */
            
            CSHIFT(VAR[mol].F[DISP][XDIR],VAR[comp].F[DISP][XDIR],
                   VAR[mol].VM[XDIR],VAR[comp].VM[XDIR],XL,BOXH,BOXL);
            CSHIFT(VAR[mol].F[DISP][YDIR],VAR[comp].F[DISP][YDIR],
                   VAR[mol].VM[YDIR],VAR[comp].VM[YDIR],YL,BOXH,BOXL);
            CSHIFT(VAR[mol].F[DISP][ZDIR],VAR[comp].F[DISP][ZDIR],
                   VAR[mol].VM[ZDIR],VAR[comp].VM[ZDIR],ZL,BOXH,BOXL);
            
            KC=0;
            for (K = 0; K < 9; K++) {
                RS[K]=XL[K]*XL[K]+YL[K]*YL[K]+ZL[K]*ZL[K];
                if (RS[K] > CUT2) 
                    KC++;
            } /* for K */
            
            if (KC != 9) {
                for (K = 0; K < 14; K++) 
                    FF[K]=0.0;
                if (RS[0] < CUT2) {
                    FF[0]=QQ4/(RS[0]*sqrt(RS[0]))+REF4;
                    LVIR = LVIR + FF[0]*RS[0];
                } /* if */
                for (K = 1; K < 5; K++) {
                    if (RS[K] < CUT2) { 
                        FF[K]= -QQ2/(RS[K]*sqrt(RS[K]))-REF2;
                        LVIR = LVIR + FF[K]*RS[K];
                    } /* if */
                    if (RS[K+4] <= CUT2) { 
                        RL[K+4]=sqrt(RS[K+4]);
                        FF[K+4]=QQ/(RS[K+4]*RL[K+4])+REF1;
                        LVIR = LVIR + FF[K+4]*RS[K+4];
                    } /* if */
                } /* for K */
                if (KC == 0) {
                    RS[9]=XL[9]*XL[9]+YL[9]*YL[9]+ZL[9]*ZL[9];
                    RL[9]=sqrt(RS[9]);
                    FF[9]=AB1*exp(-B1*RL[9])/RL[9];
                    LVIR = LVIR + FF[9]*RS[9];
                    for (K = 10; K < 14; K++) { 
                        FTEMP=AB2*exp(-B2*RL[K-5])/RL[K-5];
                        FF[K-5]=FF[K-5]+FTEMP;
                        LVIR= LVIR+FTEMP*RS[K-5];
                        RS[K]=XL[K]*XL[K]+YL[K]*YL[K]+ZL[K]*ZL[K];
                        RL[K]=sqrt(RS[K]);
                        FF[K]=(AB3*exp(-B3*RL[K])-AB4*exp(-B4*RL[K]))/RL[K];
                        LVIR = LVIR + FF[K]*RS[K];
                    } /* for K */
                } /* if KC == 0 */
                
                UPDATE_FORCES(mol, comp, XL, YL, ZL, FF, ProcID);
                
            }  /* if KC != 9 */  
        } /* for comp */
    } /* for mol */
    
    /*  accumulate the running sum from private 
        per-interaction partial sums   */
    LOCK(gl->InterfVirLock);
    *VIR = *VIR + LVIR;
    UNLOCK(gl->InterfVirLock);
    
    /* at the end of the above force-computation, comp_last */
    /* contains the number of the last molecule (no modulo) */
    /* that this process touched                            */
    
    if (comp_last > NMOL1) {
        for (mol = StartMol[ProcID]; mol < NMOL; mol++) {
            ALOCK(gl->MolLock, mol % MAXLCKS); 
            for ( dir = XDIR; dir  <= ZDIR; dir++) {
                temp_p = VAR[mol].F[DEST][dir];
                temp_p[H1] += PFORCES[ProcID][mol][dir][H1];
                temp_p[O]  += PFORCES[ProcID][mol][dir][O];
                temp_p[H2] += PFORCES[ProcID][mol][dir][H2];
            }
            AULOCK(gl->MolLock, mol % MAXLCKS);
        }
        comp = comp_last % NMOL;
        for (mol = 0; ((mol <= comp) && (mol < StartMol[ProcID])); mol++) {
            ALOCK(gl->MolLock, mol % MAXLCKS); 
            for ( dir = XDIR; dir  <= ZDIR; dir++) {
                temp_p = VAR[mol].F[DEST][dir];
                temp_p[H1] += PFORCES[ProcID][mol][dir][H1];
                temp_p[O]  += PFORCES[ProcID][mol][dir][O];
                temp_p[H2] += PFORCES[ProcID][mol][dir][H2];
            }
            AULOCK(gl->MolLock, mol % MAXLCKS);
        }
    }
    else{
        for (mol = StartMol[ProcID]; mol <= comp_last; mol++) {
            ALOCK(gl->MolLock, mol % MAXLCKS); 
            for ( dir = XDIR; dir  <= ZDIR; dir++) {
                temp_p = VAR[mol].F[DEST][dir];
                temp_p[H1] += PFORCES[ProcID][mol][dir][H1];
                temp_p[O]  += PFORCES[ProcID][mol][dir][O];
                temp_p[H2] += PFORCES[ProcID][mol][dir][H2];
            }
            AULOCK(gl->MolLock, mol % MAXLCKS);
        }
    }
    
    /* wait till all forces are updated */
    
    BARRIER(gl->InterfBar, NumProcs);
    
    /* divide final forces by masses */
    
    for (mol = StartMol[ProcID]; mol < StartMol[ProcID+1]; mol++) {
        for ( dir = XDIR; dir  <= ZDIR; dir++) {
            temp_p = VAR[mol].F[DEST][dir];
            temp_p[H1] = temp_p[H1] * FHM;
            temp_p[O]  = temp_p[O] * FOM;
            temp_p[H2] = temp_p[H2] * FHM;
        } /* for dir */
    } /* for mol */
    
}/* end of subroutine INTERF */

UPDATE_FORCES(mol, comp, XL, YL, ZL, FF, ProcID)
  /* from the computed distances etc., compute the 
     intermolecular forces and update the force (or 
     acceleration) locations */
  int mol, comp;
  double XL[], YL[], ZL[], FF[];
  unsigned ProcID;
{
    int K;
    double G110[3], G23[3], G45[3], TT1[3], TT[3], TT2[3];
    double GG[15][3];
    double *tx_p, *ty_p, *tz_p;
    
    /*   CALCULATE X-COMPONENT FORCES */
    for (K = 0; K < 14; K++)  {
        GG[K+1][XDIR] = FF[K]*XL[K];
        GG[K+1][YDIR] = FF[K]*YL[K];
        GG[K+1][ZDIR] = FF[K]*ZL[K];
    }
    
    G110[XDIR] = GG[10][XDIR]+GG[1][XDIR]*C1;
    G110[YDIR] = GG[10][YDIR]+GG[1][YDIR]*C1;
    G110[ZDIR] = GG[10][ZDIR]+GG[1][ZDIR]*C1;
    G23[XDIR] = GG[2][XDIR]+GG[3][XDIR];
    G23[YDIR] = GG[2][YDIR]+GG[3][YDIR];
    G23[ZDIR] = GG[2][ZDIR]+GG[3][ZDIR];
    G45[XDIR]=GG[4][XDIR]+GG[5][XDIR];
    G45[YDIR]=GG[4][YDIR]+GG[5][YDIR];
    G45[ZDIR]=GG[4][ZDIR]+GG[5][ZDIR];
    TT1[XDIR] =GG[1][XDIR]*C2;
    TT1[YDIR] =GG[1][YDIR]*C2;
    TT1[ZDIR] =GG[1][ZDIR]*C2;
    TT[XDIR] =G23[XDIR]*C2+TT1[XDIR];
    TT[YDIR] =G23[YDIR]*C2+TT1[YDIR];
    TT[ZDIR] =G23[ZDIR]*C2+TT1[ZDIR];
    TT2[XDIR]=G45[XDIR]*C2+TT1[XDIR];
    TT2[YDIR]=G45[YDIR]*C2+TT1[YDIR];
    TT2[ZDIR]=G45[ZDIR]*C2+TT1[ZDIR];
    /* lock locations for the molecule to be updated */
    tx_p = PFORCES[ProcID][mol][XDIR];
    ty_p = PFORCES[ProcID][mol][YDIR];
    tz_p = PFORCES[ProcID][mol][ZDIR];
    tx_p[H1] += 
        GG[6][XDIR]+GG[7][XDIR]+GG[13][XDIR]+TT[XDIR]+GG[4][XDIR];
    tx_p[O] +=
        G110[XDIR] + GG[11][XDIR] +GG[12][XDIR]+C1*G23[XDIR];
    tx_p[H2] +=
        GG[8][XDIR]+GG[9][XDIR]+GG[14][XDIR]+TT[XDIR]+GG[5][XDIR];
    ty_p[H1] += 
        GG[6][YDIR]+GG[7][YDIR]+GG[13][YDIR]+TT[YDIR]+GG[4][YDIR];
    ty_p[O]  += 
        G110[YDIR]+GG[11][YDIR]+GG[12][YDIR]+C1*G23[YDIR];
    ty_p[H2] += 
        GG[8][YDIR]+GG[9][YDIR]+GG[14][YDIR]+TT[YDIR]+GG[5][YDIR];
    tz_p[H1] +=
        GG[6][ZDIR]+GG[7][ZDIR]+GG[13][ZDIR]+TT[ZDIR]+GG[4][ZDIR];
    tz_p[O]  +=
        G110[ZDIR]+GG[11][ZDIR]+GG[12][ZDIR]+C1*G23[ZDIR];
    tz_p[H2] +=
        GG[8][ZDIR]+GG[9][ZDIR]+GG[14][ZDIR]+TT[ZDIR]+GG[5][ZDIR];
    
    tx_p = PFORCES[ProcID][comp][XDIR];
    ty_p = PFORCES[ProcID][comp][YDIR];
    tz_p = PFORCES[ProcID][comp][ZDIR];
    tx_p[H1] +=
        -GG[6][XDIR]-GG[8][XDIR]-GG[11][XDIR]-TT2[XDIR]-GG[2][XDIR];
    tx_p[O] +=
        -G110[XDIR]-GG[13][XDIR]-GG[14][XDIR]-C1*G45[XDIR];
    tx_p[H2] +=
        -GG[7][XDIR]-GG[9][XDIR]-GG[12][XDIR]-TT2[XDIR]-GG[3][XDIR];
    ty_p[H1] += 
        -GG[6][YDIR]-GG[8][YDIR]-GG[11][YDIR]-TT2[YDIR]-GG[2][YDIR];
    ty_p[O] += 
        -G110[YDIR]-GG[13][YDIR]-GG[14][YDIR]-C1*G45[YDIR];
    ty_p[H2] += 
        -GG[7][YDIR]-GG[9][YDIR]-GG[12][YDIR]-TT2[YDIR]-GG[3][YDIR];
    tz_p[H1] +=
        -GG[6][ZDIR]-GG[8][ZDIR]-GG[11][ZDIR]-TT2[ZDIR]-GG[2][ZDIR];
    tz_p[O] +=
        -G110[ZDIR]-GG[13][ZDIR]-GG[14][ZDIR]-C1*G45[ZDIR];
    tz_p[H2] +=
        -GG[7][ZDIR]-GG[9][ZDIR]-GG[12][ZDIR]-TT2[ZDIR]-GG[3][ZDIR];
}           /* end of subroutine UPDATE_FORCES */