ComputePairForce.c File Reference

#include "ComputePairForce.h"
#include <stdio.h>
#include <math.h>
#include "global.h"

Include dependency graph for ComputePairForce.c:

Go to the source code of this file.

Functions
void	ComputePairForce (int normFlag)

Function Documentation

ComputePairForce()

void ComputePairForce

(

int

normFlag

)

Definition at line 27 of file ComputePairForce.c.

                                   {
double dr[NDIM+1], fcVal, rr, ri, r, uVal;
int n, i,j;
uVal = 0.0;  uSumPair = 0.0 ; 
virSumPair = 0.0; virSumPairxx = 0.0; virSumPairyy = 0.0; virSumPairxy = 0.0;
 
for(n = 1 ; n <= nAtom ; n ++){
 fx[n] = 0.0;
 fy[n] = 0.0;
 discDragx[n] = 0.0;
 discDragy[n] = 0.0;
}
for(n = 1; n <= nPairTotal; n ++){
 PairID[n] = 0;
 Pairatom1[n] = 0;
 Pairatom2[n] = 0;
 PairXij[n] = 0.0; 
 PairYij[n] = 0.0; 
}
 
double vr[NDIM+1], fdVal, rri;
nPairActive = 0;
double meff;
meff = 0.0;
int atomIDi, atomIDj;
double atomiMass, atomjMass;
 
//int processThisPair = 1;
 
for(i=1;i<=nAtomInterface;i++){
 for(j=i+1;j<=nAtomInterface;j++){
  atomIDi = atomIDInterface[i];
  atomIDj = atomIDInterface[j];
  
  if (isBonded[atomIDi][atomIDj] == 0) { //To have pair interaction between nonbonded atoms only
  rr = 0.0; rri = 0.0; fcVal = 0.0;  fdVal = 0.0; strech = 0.0;
  RadiusIJ = 0.0;
  
  dr[1] = rx[atomIDi] - rx[atomIDj];
  if(dr[1] >= regionH[1]) 
   dr[1] -= region[1];
  else if(dr[1] < -regionH[1]) 
    dr[1] += region[1];
 
  dr[2] = ry[atomIDi] - ry[atomIDj];
   if(dr[2] >= regionH[2]){
    dr[1] -= shearDisplacement;
    if(dr[1] < -regionH[1]) dr[1] += region[1];
    dr[2] -= region[2];
  }else if(dr[2] < -regionH[2]){
   dr[1] += shearDisplacement;
    if(dr[1] >= regionH[1]) dr[1] -= region[1];
   dr[2] += region[2];
  }
 
  rr = Sqr(dr[1]) + Sqr(dr[2]);
  RadiusIJ = atomRadius[atomIDi] + atomRadius[atomIDj];
  SqrRadiusIJ = Sqr(RadiusIJ);
  if(rr < SqrRadiusIJ){
   r = sqrt(rr);
   ri = 1.0/r;
   rri = 1.0/rr;
   RadiusIJInv = 1.0/RadiusIJ;
   strech = (RadiusIJ - r);
   uVal = 0.5 * Kn * Sqr(strech);
   
   //NormFlag
   if(normFlag == 1){ 
    strech = strech * RadiusIJInv;
    uVal = 0.5 * Kn * RadiusIJ * Sqr(strech);
    }
   
   fcVal = Kn * strech * ri;
   vr[1] = vx[atomIDi] - vx[atomIDj];
   vr[2] = vy[atomIDi] - vy[atomIDj];
 
   nPairActive++;
   PairID[nPairActive] = nPairActive;
   Pairatom1[nPairActive] = atomIDi;
   Pairatom2[nPairActive] = atomIDj;
   PairXij[nPairActive] = dr[1]; 
   PairYij[nPairActive] = dr[2]; 
   
   //DampFlag = 1
   if(DampFlag == 1){
   atomiMass = atomMass[atomIDi];
   atomjMass = atomMass[atomIDj];
   meff = (atomiMass * atomjMass)/(atomiMass + atomjMass);
   fdVal = -gamman * meff * (vr[1]*dr[1] + vr[2]*dr[2]) * rri; //disc-disc drag
 
   discDragx[atomIDi] =  fdVal * dr[1]; //disc-disc drag
   discDragy[atomIDi] =  fdVal * dr[2]; //disc-disc drag
   discDragx[atomIDj] = -fdVal * dr[1]; //disc-disc drag
   discDragy[atomIDj] = -fdVal * dr[2]; //disc-disc drag
 
   discDragx[nPairActive] = discDragx[atomIDi];
   discDragy[nPairActive] = discDragy[atomIDi];
   
 
   fx[atomIDi] +=  (fcVal + fdVal) * dr[1];
   fy[atomIDi] +=  (fcVal + fdVal) * dr[2];
   fx[atomIDj] += -(fcVal + fdVal) * dr[1];
   fy[atomIDj] += -(fcVal + fdVal) * dr[2];
  }
 
  //DampFlag = 2
  else if(DampFlag == 2){
   discDragx[atomIDi] =  -gamman * vr[1]; //disc-disc drag
   discDragy[atomIDi] =  -gamman * vr[2]; //disc-disc drag
   discDragx[atomIDj] = -(-gamman * vr[1]); //disc-disc drag
   discDragy[atomIDj] = -(-gamman * vr[2]); //disc-disc drag
 
   discDragx[nPairActive] = discDragx[atomIDi];
   discDragy[nPairActive] = discDragy[atomIDi];
   
   fx[atomIDi] +=  (fcVal * dr[1] - gamman * vr[1]);
   fy[atomIDi] +=  (fcVal * dr[2] - gamman * vr[2]);
   fx[atomIDj] += -(fcVal * dr[1] - gamman * vr[1]);
   fy[atomIDj] += -(fcVal * dr[2] - gamman * vr[2]);
  }
 
  //In the following, for stress/virial term (fcVal + fdVal) is used since the total pair force = Hookean Interaction + relative velocity drag
  uSumPair +=  0.5 * uVal;
  virSumPair += (fcVal + fdVal) * rr; 
  virSumPairxx += (fcVal + fdVal) * dr[1] * dr[1];
  virSumPairyy += (fcVal + fdVal) * dr[2] * dr[2];
  virSumPairxy += (fcVal + fdVal) * dr[1] * dr[2];
    }
   }
  }
 }
}

References atomIDInterface, atomMass, atomRadius, DampFlag, discDragx, discDragy, fx, fy, gamman, isBonded, Kn, nAtom, nAtomInterface, NDIM, nPairActive, nPairTotal, Pairatom1, Pairatom2, PairID, PairXij, PairYij, RadiusIJ, RadiusIJInv, region, regionH, rx, ry, shearDisplacement, Sqr, SqrRadiusIJ, strech, uSumPair, virSumPair, virSumPairxx, virSumPairxy, virSumPairyy, vx, and vy.

Referenced by main().

Here is the caller graph for this function: