/*************************************************************************/
/*                                                                       */
/*  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 "mdvar.h"
#include "frcnst.h"
#include "water.h"
#include "wwpot.h"
#include <math.h>
#include "parameters.h"
#include "mddata.h"
#include "split.h"
#include "global.h"

POTENG(POTA,POTR,PTRF,ProcID)
  double *POTA, *POTR, *PTRF;     /* some shared sums computed by POTENG */
  unsigned ProcID;
  
{
    
    /*
      this routine calculates the potential energy of the system.  it is only
      called at periods specified by the user, typically in those
      time-steps when the user wants to print output.
      FC11 ,FC12, FC13, and FC33 are the quardratic force constants
      */
    
    int KC, K;
    int XBOX, YBOX, ZBOX;
    int X_NUM, Y_NUM, Z_NUM;
    int i, j, k;
    double R1, R2, RX, COS, DT, DR1, DR2, DR1S, DR2S, DRP, DRS;
    double XL[15], YL[15], ZL[15], RS[15], RL[15];
    double DTS;
    double LPOTA, LPOTR, LPTRF;
    struct link *curr_ptr, *neighbor_ptr;
    struct list_of_boxes *curr_box;
    double *tx_p, *ty_p, *tz_p;
    double tempa, tempb, tempc;
    
    /*  compute intra-molecular potential energy */
    
    LPOTA=0.0;
    curr_box = my_boxes[ProcID];
    while (curr_box) {
        
        i = curr_box->coord[XDIR];  /* X coordinate of box */
        j = curr_box->coord[YDIR];  /* Y coordinate of box */
        k = curr_box->coord[ZDIR];  /* Z coordinate of box */
        
        /* Go through the molecules in a box */
        
        curr_ptr = BOX[i][j][k].list;
        while (curr_ptr) {
            
            tx_p = curr_ptr->mol.F[DISP][XDIR];
            ty_p = curr_ptr->mol.F[DISP][YDIR];
            tz_p = curr_ptr->mol.F[DISP][ZDIR];
            
            curr_ptr->mol.VM[XDIR] = C1 * tx_p[O] +
                C2 * (tx_p[H1] + 
                      tx_p[H2]);
            curr_ptr->mol.VM[YDIR] = C1*ty_p[O] +
                C2*(ty_p[H1] + 
                    ty_p[H2]);
            curr_ptr->mol.VM[ZDIR] = C1*tz_p[O] +
                C2*(tz_p[H1] + 
                    tz_p[H2]);
            tempa = tx_p[O] - tx_p[H1];
            tempb = ty_p[O] - ty_p[H1];
            tempc = tz_p[O] - tz_p[H1];
            R1 = tempa * tempa + tempb * tempb + tempc * tempc;
            tempa = tx_p[O] - tx_p[H2];
            tempb = ty_p[O] - ty_p[H2];
            tempc = tz_p[O] - tz_p[H2];
            R2 = tempa * tempa + tempb * tempb + tempc * tempc;
            
            RX = ((tx_p[O] - tx_p[H1]) *
                  (tx_p[O] - tx_p[H2])) +
                      ((ty_p[O] - ty_p[H1]) *
                       (ty_p[O] - ty_p[H2])) +
                           ((tz_p[O] - tz_p[H1]) *
                            (tz_p[O] - tz_p[H2]));
            
            R1=sqrt(R1);
            R2=sqrt(R2);
            COS=RX/(R1*R2);
            DT=(acos(COS)-ANGLE)*ROH;
            DR1=R1-ROH;
            DR2=R2-ROH;
            DR1S=DR1*DR1;
            DR2S=DR2*DR2;
            DRP=DR1+DR2;
            DTS=DT*DT;
            LPOTA += (FC11*(DR1S+DR2S)+FC33*DTS)*0.5
                +FC12*DR1*DR2+FC13*DRP*DT
                    +(FC111*(DR1S*DR1+DR2S*DR2)+FC333*DTS*DT+FC112*DRP*DR1*DR2+
                      FC113*(DR1S+DR2S)*DT+FC123*DR1*DR2*DT+FC133*DRP*DTS)*ROHI;
            
            LPOTA += (FC1111*(DR1S*DR1S+DR2S*DR2S)+FC3333*DTS*DTS+
                      FC1112*(DR1S+DR2S)*DR1*DR2+FC1122*DR1S*DR2S+
                      FC1113*(DR1S*DR1+DR2S*DR2)*DT+FC1123*DRP*DR1*DR2*DT+
                      FC1133*(DR1S+DR2S)*DTS+FC1233*DR1*DR2*DTS+
                      FC1333*DRP*DTS*DT)*ROHI2;
            
            curr_ptr = curr_ptr->next_mol;
        } /* while curr_ptr */
        curr_box = curr_box->next_box;
    } /* while curr_box */
    
    
    BARRIER(gl->PotengBar, NumProcs);
    
    /*  compute inter-molecular potential energy */
    
    LPOTR=0.0;
    LPTRF=0.0;
    curr_box = my_boxes[ProcID];
    while (curr_box) {
        
        i = curr_box->coord[XDIR];  /* X coordinate of box */
        j = curr_box->coord[YDIR];  /* Y coordinate of box */
        k = curr_box->coord[ZDIR];  /* Z coordinate of box */
        
        /* loop over nearest neighbor boxes */
        
        for (XBOX=i-1; XBOX<=i+1; XBOX++) {
            for (YBOX=j-1; YBOX<=j+1; YBOX++) {
                for (ZBOX=k-1; ZBOX<=k+1; ZBOX++) {
                    
                    X_NUM = XBOX;
                    Y_NUM = YBOX;
                    Z_NUM = ZBOX;
                    
                    if ((BOX_PER_SIDE == 2) && ((XBOX < 0) || (XBOX == 2) ||
                                                (YBOX < 0) || (YBOX == 2) || 
                                                (ZBOX < 0) || (ZBOX == 2)))
                        continue;
                    
                    /* Make box number valid if out of box */
                    
                    if (X_NUM == -1)
                        X_NUM += BOX_PER_SIDE;
                    else if (X_NUM >= BOX_PER_SIDE)
                        X_NUM -= BOX_PER_SIDE;
                    if (Y_NUM == -1)
                        Y_NUM += BOX_PER_SIDE;
                    else if (Y_NUM >= BOX_PER_SIDE)
                        Y_NUM -= BOX_PER_SIDE;
                    if (Z_NUM == -1)
                        Z_NUM += BOX_PER_SIDE;
                    else if (Z_NUM >= BOX_PER_SIDE)
                        Z_NUM -= BOX_PER_SIDE;
                    
                    
                    /* Don't do current box more than once */
                    if ((X_NUM == i) && (Y_NUM == j) && (Z_NUM == k) &&
                        ((XBOX != i) || (YBOX != j) || (ZBOX !=k))) {
                        continue;
                    }
                    
                    neighbor_ptr = BOX[X_NUM][Y_NUM][Z_NUM].list;
                    while (neighbor_ptr) {
                        /* go through current box list */
                        curr_ptr = BOX[i][j][k].list;
                        while (curr_ptr) {
                            
                            /* Don't do interaction with same molecule */
                            
                            if (curr_ptr == neighbor_ptr) {
                                curr_ptr = curr_ptr->next_mol;
                                continue;
                            }
                            
                            /*  compute some intermolecular distances */
                            CSHIFT(curr_ptr->mol.F[DISP][XDIR],neighbor_ptr->mol.F[DISP][XDIR],
                                   curr_ptr->mol.VM[XDIR],neighbor_ptr->mol.VM[XDIR],XL,BOXH,BOXL);
                            CSHIFT(curr_ptr->mol.F[DISP][YDIR],neighbor_ptr->mol.F[DISP][YDIR],
                                   curr_ptr->mol.VM[YDIR],neighbor_ptr->mol.VM[YDIR],YL,BOXH,BOXL);
                            CSHIFT(curr_ptr->mol.F[DISP][ZDIR],neighbor_ptr->mol.F[DISP][ZDIR],
                                   curr_ptr->mol.VM[ZDIR],neighbor_ptr->mol.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 < 9; K++) {
                                    if (RS[K] <= CUT2) {
                                        RL[K]=sqrt(RS[K]);
                                    }
                                    else {
                                        RL[K]=CUTOFF;
                                        RS[K]=CUT2;
                                    } /* else */
                                } /* for K */
                                
                                LPOTR= LPOTR-QQ2/RL[1]-QQ2/RL[2]-QQ2/RL[3]-QQ2/RL[4]
                                    + QQ/RL[5]+ QQ/RL[6]+ QQ/RL[7]+ QQ/RL[8]
                                        + QQ4/RL[0];
                                
                                LPTRF= LPTRF-REF2*RS[0]-REF1*((RS[5]+RS[6]+RS[7]+RS[8])*0.5
                                                              -RS[1]-RS[2]-RS[3]-RS[4]);
                                
                                if (KC <= 0) {
                                    for (K = 9; K <  14; K++)  {
                                        RL[K]=sqrt(XL[K]*XL[K]+YL[K]*YL[K]+ZL[K]*ZL[K]);
                                    }
                                    
                                    LPOTR= LPOTR+A1* exp(-B1*RL[9])
                                        +A2*(exp(-B2*RL[ 5])+exp(-B2*RL[ 6])
                                             +exp(-B2*RL[ 7])+exp(-B2*RL[ 8]))
                                            +A3*(exp(-B3*RL[10])+exp(-B3*RL[11])
                                                 +exp(-B3*RL[12])+exp(-B3*RL[13]))
                                                -A4*(exp(-B4*RL[10])+exp(-B4*RL[11])
                                                     +exp(-B4*RL[12])+exp(-B4*RL[13]));
                                } /* if KC <= 0 */
                                
                                
                            } /* if KC != 9 */
                            curr_ptr = curr_ptr->next_mol;
                        }
                        neighbor_ptr = neighbor_ptr->next_mol;
                    }
                }
            }
        } /* neighbor boxes for loops */
        curr_box = curr_box->next_box;
    } /* while curr_box */
    
    LPOTR = LPOTR/2.0;
    LPTRF = LPTRF/2.0;
    
    /* update shared sums from computed private sums */
    
    LOCK(gl->PotengSumLock);
    *POTA = *POTA + LPOTA;
    *POTR = *POTR + LPOTR;
    *PTRF = *PTRF + LPTRF;
    UNLOCK(gl->PotengSumLock);
    
} /* end of subroutine POTENG */