main.c File Reference

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdbool.h>
#include <time.h>
#include "global.h"
#include "ComputeBondForce.h"
#include "ComputePairForce.h"

Include dependency graph for main.c:

Go to the source code of this file.

Macros
#define	DEFINE_GLOBALS

Functions
void	Init ()
void	SetupJob ()
void	EvalSpacetimeCorr ()
void	Trajectory ()
void	DumpState ()
void	ComputeForcesCells ()
void	ApplyBoundaryCond ()
void	EvalProps ()
void	AccumProps (int icode)
void	PrintSummary ()
void	PrintVrms ()
void	VelocityVerletStep (int icode)
void	ApplyForce ()
void	ApplyLeesEdwardsBoundaryCond ()
void	PrintStress ()
void	Close ()
void	PrintMomentum ()
void	DisplaceAtoms ()
void	DumpRestart ()
bool	HaltConditionCheck (double value)
void	EvalCom ()
void	PrintCom ()
void	EvalVrms ()
void	EvalUnwrap ()
void	DumpBonds ()
void	DumpPairs ()
void	WriteBinaryRestart ()
void	PrintForceSum ()
int	main (int argc, char **argv)

Variables
char *	prefix = NULL

Macro Definition Documentation

DEFINE_GLOBALS

#define DEFINE_GLOBALS

Definition at line 7 of file main.c.

Function Documentation

AccumProps()

void AccumProps

(

int

icode

)

Definition at line 25 of file AccumProps.c.

                          {
 if(icode == 0){
 sPotEnergy = ssPotEnergy = 0.;
 sKinEnergy = ssKinEnergy = 0.;
 sPressure = ssPressure = 0.;
 sTotEnergy = ssTotEnergy = 0.; 
 svirSum = 0.;
 }else if(icode == 1){
 sPotEnergy += potEnergy;
 ssPotEnergy += Sqr(potEnergy);
 sKinEnergy += kinEnergy;
 ssKinEnergy += Sqr(kinEnergy);
 sTotEnergy += totEnergy;
 ssTotEnergy += Sqr(totEnergy);
 sPressure += pressure;
 ssPressure += Sqr(pressure);
 svirSum += virSum;
 }else if(icode == 2){
 sPotEnergy /= stepAvg;
 ssPotEnergy /= sqrt(ssPotEnergy/stepAvg - Sqr(sPotEnergy));
 sTotEnergy /= stepAvg;
 ssTotEnergy = sqrt(ssTotEnergy/stepAvg - Sqr(sTotEnergy));
 sKinEnergy /= stepAvg;
 ssKinEnergy = sqrt(ssKinEnergy/stepAvg - Sqr(sKinEnergy));
 sPressure /= stepAvg;
 ssPressure = sqrt(ssPressure/stepAvg - Sqr(sPressure));
 svirSum /= stepAvg;
 } }

ApplyBoundaryCond()

void ApplyBoundaryCond

(

)

Definition at line 27 of file ApplyBoundaryCond.c.

                        {
  int n;
  for(n = 1 ; n <= nAtom ; n ++){
   if(strcmp(xBoundary, "p") == 0 && strcmp(yBoundary, "p") == 0){         // P.B.C along x and y axis
    rx[n] -= region[1]*rint(rx[n]/region[1]);
    ry[n] -= region[2]*rint(ry[n]/region[2]);
   } else if (strcmp(xBoundary, "r") == 0 && strcmp(yBoundary, "r") == 0){   //R.B.C. along x and y axis
     if((rx[n] + atomRadius[n]) >= regionH[1]){
        rx[n] = 0.999999*regionH[1] - atomRadius[n]; vx[n] = -vx[n] ;
     }if((rx[n]-atomRadius[n]) < -regionH[1]){
        rx[n] = -0.999999*regionH[1] + atomRadius[n]; vx[n] = -vx[n] ;
    }
     if((ry[n] + atomRadius[n])>= regionH[2]){
        ry[n] = 0.999999*regionH[2] - atomRadius[n]; vy[n] = -vy[n] ;
     }if((ry[n]-atomRadius[n]) < -regionH[2]){
        ry[n] = -0.999999*regionH[2] + atomRadius[n]; vy[n] = -vy[n] ;
    }} 
    else if (strcmp(xBoundary, "p") == 0 && strcmp(yBoundary, "r") == 0){   //P.B.C. along x and R.B.C along y axis
    rx[n] -= region[1]*rint(rx[n]/region[1]); 
    if((ry[n] + atomRadius[n]) >= regionH[2]){
      ry[n] = 0.999999*regionH[2] - atomRadius[n]; vy[n] = -vy[n] ;
     }if((ry[n] - atomRadius[n]) < -regionH[2]){
        ry[n] = -0.999999*regionH[2] + atomRadius[n]; vy[n] = -vy[n] ;
  }} 
    else if(strcmp(xBoundary, "r") == 0 && strcmp(yBoundary, "p") == 0){   //R.B.C. along x and P.B.C along y axis
    if((rx[n] + atomRadius[n]) >= regionH[1]){
       rx[n] = 0.999999*regionH[1] - atomRadius[n]; vx[n] = -vx[n] ;
     }if((rx[n] - atomRadius[n]) < -regionH[1]){
      rx[n] = -0.999999*regionH[1] + atomRadius[n]; vx[n] = -vx[n] ;
     } 
     ry[n] -= region[2]*rint(ry[n]/region[2]);  
  } else {
    // Print error message and exit the program
    fprintf(fpresult, "Error: Invalid boundary configuration: '%s %s'\n", xBoundary, yBoundary);
    exit(EXIT_FAILURE); // Exit with failure status 
  }
 }
}

References atomRadius, fpresult, nAtom, region, regionH, rx, ry, vx, vy, xBoundary, and yBoundary.

Referenced by main().

Here is the caller graph for this function:

ApplyForce()

void ApplyForce

(

)

Definition at line 25 of file ApplyForce.c.

                 {
 int n;
 double lx;
 lx = regionH[1];  
 fyExtern =  (FyBylx * lx)/nAtomBlock;
 fxExtern =  fxByfy * fyExtern;
 forceSumxExtern = fxExtern*nAtomBlock;  forceSumyExtern = fyExtern*nAtomBlock;
 
 for(n = 1; n <= nAtom; n ++){
 if(molID[n] == 2){
  fx[n] += fxExtern; 
  fy[n] -= fyExtern;
} } }

References forceSumxExtern, forceSumyExtern, fx, fxByfy, fxExtern, fy, FyBylx, fyExtern, molID, nAtom, nAtomBlock, and regionH.

Referenced by main().

Here is the caller graph for this function:

ApplyLeesEdwardsBoundaryCond()

void ApplyLeesEdwardsBoundaryCond

(

)

Definition at line 25 of file ApplyLeesEdwardsBoundaryCond.c.

                                    {
 int n;
 for (n = 1; n <= nAtom; n++) {
//PBC along x-direction
 if(rx[n] >= regionH[1])
   rx[n] -= region[1];
 else if(rx[n] < -regionH[1])
   rx[n] += region[1];
   
//LEBC along y-direction
  if(ry[n] >= regionH[2]){
   rx[n] -= shearDisplacement;
   if(rx[n] < -regionH[1]) rx[n] += region[1];
   //vx[n] -= shearVelocity; 
   ry[n] -= region[2];
  }else if(ry[n] < -regionH[2]){
   rx[n] += shearDisplacement;
   if(rx[n] >= regionH[1]) rx[n] -= region[1];
   //vx[n] += shearVelocity; 
   ry[n] += region[2];
  }
 } 
}

References nAtom, region, regionH, rx, ry, and shearDisplacement.

Close()

void Close

(

)

Definition at line 25 of file Close.c.

            {
  int n;
  free(rx); free(ry); free(vx); free(vy); free(ax); free(ay); free(fx); free(fy);
  free(fax);
  free(fay);
  free(cellList);
 
  free(atomID); free(atomType); free(atomRadius); free(atomMass);
  free(speed);
  free(atom1); free(atom2); free(BondID);
  free(BondType); free(kb); free(ro);
  free(ImageX); free(ImageY); free(rxUnwrap); free(ryUnwrap);
  free(atomIDInterface);
  free(PairID); free(Pairatom1); free(Pairatom2); 
  free(PairXij); free(PairYij);
  free(molID);
  
  for (n = 0; n <= nAtom; n++) {
   free(isBonded[n]);
   }
  free(isBonded);
 
  for (n = 0; n <= nBuffCorr; n++){
  free(cfOrg[n]);
  free(spacetimeCorr[n]);
  }
  free(cfOrg);
  free(spacetimeCorr);
  free(cfVal);
  free(indexCorr);
  free(spacetimeCorrAv);
 
  free(indexAcf);
  free(viscAcfOrg);
  free(viscAcfAv);
  for(n = 0 ; n <= nBuffAcf ; n ++)
   free(viscAcf[n]);
  free(viscAcf);
 
 }

References atom1, atom2, atomID, atomIDInterface, atomMass, atomRadius, atomType, ax, ay, BondID, BondType, cellList, cfOrg, cfVal, fax, fay, fx, fy, ImageX, ImageY, indexAcf, indexCorr, isBonded, kb, molID, nAtom, nBuffAcf, nBuffCorr, Pairatom1, Pairatom2, PairID, PairXij, PairYij, ro, rx, rxUnwrap, ry, ryUnwrap, spacetimeCorr, spacetimeCorrAv, speed, viscAcf, viscAcfAv, viscAcfOrg, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

ComputeForcesCells()

void ComputeForcesCells

(

)

Definition at line 25 of file ComputeForcesCells.c.

                         {
  double dr[NDIM+1], invWid[NDIM+1], shift[NDIM+1], f, fcVal, rr, ri, r, uVal;
  int c, I, J, m1, m1X, m1Y, m2, m2X, m2Y, n, offset;
  int iofX[] = {0, 0, 1, 1, 0, -1, -1, -1, 0, 1},
      iofY[] = {0, 0, 0, 1 ,1, 1, 0, -1, -1, -1};
 
  invWid[1] = cells[1]/region[1];
  invWid[2] = cells[2]/region[2];
 
  for(n = nAtom+1; n <= nAtom+cells[1]*cells[2] ; n++)
    cellList[n] = 0;
  
  for(n = 1 ; n <= nAtom ; n ++){
    c = ((int)((ry[n] + regionH[2])*invWid[2]))*cells[1] + (int)((rx[n]+regionH[1])*invWid[1]) + nAtom+ 1;
    cellList[n] = cellList[c];
    cellList[c] = n;
  }
  
  for(n = 1 ; n <= nAtom ; n ++){
    ax[n] = 0.;
    ay[n] = 0.;
  }
  
  uSum = 0.0 ;
  virSum = 0.0;
  rfAtom = 0.0;
  RadiusIJ = 0.0;
 
  gamman = 1.0;
  double vr[NDIM+1], fd, fdVal, rrinv;
  rrinv = 0.0;
  fd = 0.0;
  fdVal = 0.0; 
  
  int start = 1 + rank*(cells[2]/size);
  int end = (rank+1)*(cells[2]/size);
 
  for(m1Y = start ; m1Y <= end ; m1Y ++){
    for(m1X = 1 ; m1X <= cells[1] ; m1X ++){
      m1 = (m1Y-1) * cells[1] + m1X + nAtom;
      for(offset = 1 ; offset <= 9 ; offset ++){
    m2X = m1X + iofX[offset]; shift[1] = 0.;
    if(m2X > cells[1]){
      m2X = 1; shift[1] = region[1];
    }else if(m2X == 0){
      m2X = cells[1]; shift[1] = -region[1];
    }
    m2Y = m1Y + iofY[offset]; shift[2] = 0.;
    if(m2Y > cells[2]){
      m2Y = 1; shift[2] = region[2];
    }else if(m2Y == 0){
      m2Y = cells[2]; shift[2] = -region[2];
    }
    m2 = (m2Y-1)*cells[1] + m2X + nAtom;
    I = cellList[m1];
    while(I > 0){
      J = cellList[m2];
      while(J > 0){
        if(m1 == m2 && J != I && (atomRadius[I] > 0. && atomRadius[J] > 0.)){
          dr[1] = rx[I] - rx[J] - shift[1];
          dr[2] = ry[I] - ry[J] - shift[2];
          rr = Sqr(dr[1]) + Sqr(dr[2]);
          RadiusIJ = atomRadius[I] + atomRadius[J];
          SqrRadiusIJ = Sqr(RadiusIJ);
          if(rr < SqrRadiusIJ){
        r = sqrt(rr);
        ri = 1.0/r;
                rrinv = 1.0/rr;
                vr[1] = vx[I] - vx[J];
                vr[2] = vy[I] - vy[J];
        RadiusIJInv = 1.0/RadiusIJ;
        uVal = Sqr(1.0 - r * RadiusIJInv);
        fcVal = 2.0 * RadiusIJInv * (1.0 - r * RadiusIJInv) *ri;
                fdVal = -gamman * (vr[1]*dr[1] + vr[2]*dr[2]) * rrinv; //disc-disc drag
 
        f  = fcVal * dr[1];
                fd = fdVal * dr[1];
        ax[I] += (f + fd);
                discDragx[I] += fd; //disc-disc drag
 
        f = fcVal * dr[2];
                fd = fdVal * dr[2];
        ay[I] += (f + fd);
                discDragy[I] += fd; //disc-disc drag
 
        uSum +=  0.5 * uVal;
        virSum += 0.5 * fcVal * rr;
        rfAtom += 0.5 * dr[1] * fcVal * dr[2];
          }
        }else if(m1 != m2 && (atomRadius[I] > 0. && atomRadius[J] > 0.)){
          dr[1] = rx[I] - rx[J] - shift[1];
          dr[2] = ry[I] - ry[J] - shift[2];
          rr = Sqr(dr[1]) + Sqr(dr[2]);
          RadiusIJ = atomRadius[I] + atomRadius[J];
          SqrRadiusIJ = Sqr(RadiusIJ);
          if(rr < SqrRadiusIJ){
        r = sqrt(rr);
        ri = 1.0/r;
                rrinv = 1.0/r;
                vr[1] = vx[I] - vx[J];
                vr[2] = vy[I] - vy[J];
        RadiusIJInv = 1.0/RadiusIJ;
        uVal = Sqr(1.0 - r * RadiusIJInv);
        fcVal = 2.0 * RadiusIJInv * (1.0 - r * RadiusIJInv) *ri;
                fdVal = -gamman * (vr[1]*dr[1] + vr[2]*dr[2]) * rrinv; //disc-disc drag
 
        f  = fcVal * dr[1];
                fd =  fdVal * dr[1];
        ax[I] += (f + fd);
                discDragx[I] += fd; //disc-disc drag
 
        f  = fcVal * dr[2];
                fd = fdVal * dr[2];
        ay[I] += (f + fd);
                discDragy[I] += fd; //disc-disc drag
 
        uSum +=  0.5 * uVal;
        virSum += 0.5 * fcVal * rr;
        rfAtom += 0.5 * dr[1] * fcVal * dr[2];
          }
        }
            J = cellList[J];
      }
      I = cellList[I];
    }
      }
    }
  }
}

References atomRadius, ax, ay, cellList, cells, discDragx, discDragy, gamman, nAtom, NDIM, RadiusIJ, RadiusIJInv, rank, region, regionH, rfAtom, rx, ry, size, Sqr, SqrRadiusIJ, uSum, virSum, vx, and vy.

DisplaceAtoms()

void DisplaceAtoms

(

)

Definition at line 25 of file DisplaceAtoms.c.

                    {
 int n;
  for(n = 1; n <= nAtom; n ++){
   if(molID[n] == 2){
    rx[n] += DeltaX; 
    ry[n] += DeltaY;
} } }

References DeltaX, DeltaY, molID, nAtom, rx, and ry.

Referenced by main().

Here is the caller graph for this function:

DumpBonds()

void DumpBonds

(

)

Definition at line 24 of file DumpBonds.c.

                {
  int n;
  //Trajectory file in LAMMPS dump format for OVITO visualization
  fprintf(fpbond, "ITEM: TIMESTEP\n");
  fprintf(fpbond, "%lf\n",timeNow);
  fprintf(fpbond, "ITEM: NUMBER OF ENTRIES\n");
  fprintf(fpbond, "%d\n",nBond);
  fprintf(fpbond, "ITEM: BOX BOUNDS pp ff pp\n");
  fprintf(fpbond, "%lf %lf xlo xhi\n", -regionH[1], regionH[1]);
  fprintf(fpbond, "%lf %lf ylo yhi\n", -regionH[2], regionH[2]);
  fprintf(fpbond, "%lf %lf zlo zhi\n", -0.1, 0.1);
  fprintf(fpbond, "ITEM: ENTRIES BondID, BondType, atom1 atom2 BondLength BondLengthEqul nodeDragx1 nodeDragy1\n");
 
  for(n=1; n<=nBond; n++)
   fprintf(fpbond, "%d %d %d %d %0.16lf %0.16lf %0.16lf %0.16lf\n", BondID[n], BondType[n], atom1[n], atom2[n], 
   BondLength[n], ro[n], nodeDragx[atom1[n]], nodeDragy[atom1[n]]);
  }

References atom1, atom2, BondID, BondLength, BondType, fpbond, nBond, nodeDragx, nodeDragy, regionH, ro, and timeNow.

Referenced by main().

Here is the caller graph for this function:

DumpPairs()

void DumpPairs

(

)

Definition at line 25 of file DumpPairs.c.

                {
  int n;
  //Trajectory file in LAMMPS dump format for OVITO visualization
  fprintf(fppair, "ITEM: TIMESTEP\n");
  fprintf(fppair, "%lf\n",timeNow);
  fprintf(fppair, "ITEM: NUMBER OF ENTRIES\n");
  fprintf(fppair, "%d\n",nPairActive);
  fprintf(fppair, "ITEM: BOX BOUNDS pp ff pp\n");
  fprintf(fppair, "%lf %lf xlo xhi\n", -regionH[1], regionH[1]);
  fprintf(fppair, "%lf %lf ylo yhi\n", -regionH[2], regionH[2]);
  fprintf(fppair, "%lf %lf zlo zhi\n", -0.1, 0.1);
  fprintf(fppair, "ITEM: ENTRIES index, atom1 atom2 xij yij discDragx1 discDragy1\n");
 
  for(n=1; n<=nPairActive; n++)
   fprintf(fppair, "%d %d %d %0.16lf %0.16lf %0.16lf %0.16lf\n", PairID[n], Pairatom1[n], Pairatom2[n], 
   PairXij[n], PairYij[n], discDragx[n], discDragy[n]);
 
  }

References discDragx, discDragy, fppair, nPairActive, Pairatom1, Pairatom2, PairID, PairXij, PairYij, regionH, and timeNow.

Referenced by main().

Here is the caller graph for this function:

DumpRestart()

void DumpRestart

(

)

Definition at line 25 of file DumpRestart.c.

                   {
 char DUMP[256];
 FILE *fpDUMP;
 sprintf(DUMP, "%s.Restart", prefix);
 fpDUMP = fopen(DUMP, "w");
 if(fpDUMP == NULL) {
  fprintf(stderr, "Error opening file %s for writing\n", DUMP);
  return;
 }
 
  fprintf(fpDUMP, "timeNow %lf\n", timeNow);
  fprintf(fpDUMP, "nAtom %d\n",  nAtom);
  fprintf(fpDUMP, "nBond %d\n", nBond);
  fprintf(fpDUMP, "nAtomType %d\n", nAtomType);
  fprintf(fpDUMP, "nBondType %d\n", nBondType);
  fprintf(fpDUMP, "region[1] %0.16lf\n", region[1]);
  fprintf(fpDUMP, "region[2] %0.16lf\n", region[2]);
 
  int n; 
  fprintf(fpDUMP, "Atoms\n");
  for(n = 1; n <= nAtom; n ++)
   fprintf(fpDUMP, "%d %d %d %0.2lf %0.16lf %0.16lf %0.16lf %0.16lf\n", atomID[n], molID[n], atomType[n], atomRadius[n], rx[n], ry[n], vx[n], vy[n]);
  
 
  fprintf(fpDUMP, "Bonds\n"); 
  for(n=1; n<=nBond; n++)
  fprintf(fpDUMP,  "%d %d %d %d %0.2lf %0.16lf\n", BondID[n], BondType[n], atom1[n], atom2[n], kb[n], ro[n]);
 
  fclose(fpDUMP);
}

References atom1, atom2, atomID, atomRadius, atomType, BondID, BondType, kb, molID, nAtom, nAtomType, nBond, nBondType, prefix, region, ro, rx, ry, timeNow, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

DumpState()

void DumpState

(

)

Definition at line 25 of file DumpState.c.

                 {
 char DUMP[256];
 FILE *fpDUMP;
 sprintf(DUMP, "%s.STATE", prefix);
 fpDUMP = fopen(DUMP, "w");
 if(fpDUMP == NULL) {
  fprintf(stderr, "Error opening file %s for writing\n", DUMP);
  return;
 }
 
  fprintf(fpDUMP, "ITEM: TIMESTEP\n");
  fprintf(fpDUMP, "%lf\n",timeNow);
  fprintf(fpDUMP, "ITEM: NUMBER OF ATOMS\n");
  fprintf(fpDUMP, "%d\n",nAtom);
  fprintf(fpDUMP, "ITEM: BOX BOUNDS pp pp pp\n");
  fprintf(fpDUMP, "%lf %lf xlo xhi\n", -regionH[1], regionH[1]);
  fprintf(fpDUMP, "%lf %lf ylo yhi\n", -regionH[2], regionH[2]);
  fprintf(fpDUMP, "%lf %lf zlo zhi\n",  -0.1, 0.1);
  fprintf(fpDUMP, "ITEM: ATOMS id mol type radius x y vx vy fx fy\n");
  int n;
  for (n = 1; n <= nAtom; n++) {
   fprintf(fpDUMP, "%d\t %d\t %d\t %0.2lf\t %0.16lf\t %0.16lf\t %0.16lf\t %0.16lf\t %0.16lf\t %0.16lf\n", 
   atomID[n], molID[n], atomType[n], atomRadius[n], rx[n], ry[n], vx[n], vy[n], ax[n], ay[n]);
  }
  fclose(fpDUMP);
}

References atomID, atomRadius, atomType, ax, ay, molID, nAtom, prefix, regionH, rx, ry, timeNow, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

EvalCom()

void EvalCom

(

)

Definition at line 27 of file EvalCom.c.

              {
 int n;
 ComX = 0.0; ComY = 0.0; ComXRatio = 0.0; ComYRatio = 0.0;
 TotalMass = 0.0;
  
 for(n=1; n<=nAtom; n++){
 if(molID[n] == 2){
  ComX += atomMass[n] * rxUnwrap[n]; 
  ComY += atomMass[n] * ryUnwrap[n]; 
  TotalMass += atomMass[n];
  } }
 
  ComX = ComX/TotalMass;
  ComY = ComY/TotalMass;
 
  if(timeNow == 0.0){
  ComX0 = ComX; ComY0 = ComY;
  }
  ComXRatio = ComX/ComX0;   ComYRatio = ComY/ComY0;
 }

References atomMass, ComX, ComX0, ComXRatio, ComY, ComY0, ComYRatio, molID, nAtom, rxUnwrap, ryUnwrap, timeNow, and TotalMass.

Referenced by main().

Here is the caller graph for this function:

EvalProps()

void EvalProps

(

)

Definition at line 26 of file EvalProps.c.

                 {
 double v;
 int n;
 double atomMassn;
 double KineEnrXSum, KineEnrYSum;
 virSum = 0.0;
 vSumX = 0.0; vSumY = 0.0; vSum = 0.0; vvSum = 0.0;
 KineEnrXSum = 0.0; KineEnrYSum = 0.0;
 
 for (n = 1; n <= nAtom; n++) {
  // Initialize v with a default value to avoid "uninitialized" warning.
  v = 0.0;
  atomMassn = atomMass[n];
  // X direction velocity
  if (strcmp(solver, "Verlet") == 0) {
    v = vx[n];
  } else if (strcmp(solver, "LeapFrog") == 0) {
    v = vx[n] - 0.5 * deltaT * ax[n];
  }
   vSum += v;
   vvSum += Sqr(v);
   KineEnrXSum += 0.5 * atomMassn * Sqr(v);
   vSumX += v; 
   // Y direction velocity
   if (strcmp(solver, "Verlet") == 0) {
   v = vy[n];
   } else if (strcmp(solver, "LeapFrog") == 0) {
   v = vy[n] - 0.5 * deltaT * ay[n];
   }
   vSum += v;
   vSumY += v;
   vvSum += Sqr(v);
   KineEnrYSum += 0.5 * atomMassn * Sqr(v);
  }
 
  kinEnergy = (KineEnrXSum + KineEnrYSum) / nAtom ;
  uSumPairPerAtom = uSumPair / nAtom ;
  BondEnergyPerAtom = TotalBondEnergy / (0.5*nAtom); //Factor of 0.5 since each atom has one half the bond energy
  potEnergy = uSumPairPerAtom +  BondEnergyPerAtom ;
  totEnergy = kinEnergy + potEnergy;
  virSumxx = virSumPairxx  + virSumBondxx ;
  virSumyy = virSumPairyy  + virSumBondyy ;
  virSumxy = virSumPairxy  + virSumBondxy ;
  virSum = virSumPair + virSumBond;
  pressure = density * (vvSum + virSum) / (nAtom * NDIM);
 
}

References atomMass, ax, ay, BondEnergyPerAtom, deltaT, density, kinEnergy, nAtom, NDIM, potEnergy, pressure, solver, Sqr, TotalBondEnergy, totEnergy, uSumPair, uSumPairPerAtom, virSum, virSumBond, virSumBondxx, virSumBondxy, virSumBondyy, virSumPair, virSumPairxx, virSumPairxy, virSumPairyy, virSumxx, virSumxy, virSumyy, vSum, vSumX, vSumY, vvSum, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

EvalSpacetimeCorr()

void EvalSpacetimeCorr

(

)

Definition at line 26 of file EvalSpacetimeCorr.c.

                         {
  real cosV=0., cosV0=0., cosV1=0., cosV2=0., sinV=0., sinV1=0., sinV2=0.;
  real COSA, SINA, COSV, SINV;
  int j, m, n, nb, ni, nv;
  real kMin = 2. * M_PI / region[1];
  real kMax = M_PI;
  real deltaK = (kMax - kMin) / nFunCorr;
  
  for (j = 1; j <= 2*nFunCorr; j++)
    cfVal[j] = 0.;
 
  for (n = 1; n <= nAtom; n++){
    j = 1;
    COSA = cos(kMin*rx[n]);
    SINA = sin(kMin*rx[n]);
    for (m = 1; m <= nFunCorr; m++){
      if(m == 1){
    cosV = cos(deltaK*rx[n]);
    sinV = sin(deltaK*rx[n]);
    cosV0 = cosV;
      }else if(m == 2){
    cosV1 = cosV;
    sinV1 = sinV;
    cosV = 2.*cosV0*cosV1-1;
    sinV = 2.*cosV0*sinV1;
      }else{
    cosV2 = cosV1;
    sinV2 = sinV1;
    cosV1 = cosV;
    sinV1 = sinV;
    cosV = 2.*cosV0*cosV1-cosV2;
    sinV = 2.*cosV0*sinV1-sinV2;
      }
      COSV = COSA*cosV - SINA*sinV;
      SINV = SINA*cosV + COSA*sinV;
      cfVal[j] += COSV;
      cfVal[j+1] += SINV;
      j += 2;
    }
  }
 
  for (nb = 1; nb <= nBuffCorr; nb++){
    indexCorr[nb] += 1;
    if (indexCorr[nb] <= 0) continue;
    ni = nFunCorr * (indexCorr[nb] - 1);
    if (indexCorr[nb] == 1){
      for (j = 1; j <= 2*nFunCorr; j++)
    cfOrg[nb][j] = cfVal[j];
    }
 
    for (j = 1; j <= nFunCorr; j++)
      spacetimeCorr[nb][ni + j] = 0.;
    
    j = 1;
    for (m = 1; m <= nFunCorr; m++){
      nv = m + ni;
      spacetimeCorr[nb][nv] += cfVal[j] * cfOrg[nb][j] + cfVal[j + 1] * cfOrg[nb][j + 1];
      j += 2;
    }
 
  }
  
  // ACCUMULATE SPACETIME CORRELATIONS
  for (nb = 1; nb <= nBuffCorr; nb++){
   if (indexCorr[nb] == nValCorr){
     for (j = 1; j <= nFunCorr*nValCorr; j++)
       spacetimeCorrAv[j] += spacetimeCorr[nb][j];
     indexCorr[nb] = 0.;
     countCorrAv ++;
     if (countCorrAv == limitCorrAv){
       for (j = 1; j <= nFunCorr*nValCorr; j++)
     spacetimeCorrAv[j] /= (nAtom*limitCorrAv);
       fprintf(fpdnsty,"NDIM %d\n", NDIM);
       fprintf(fpdnsty,"nAtom %d\n", nAtom);
       fprintf(fpdnsty,"region %lf\n", region[1]);
       fprintf(fpdnsty,"nFunCorr %d\n", nFunCorr);
       fprintf(fpdnsty,"limitCorrAv %d\n", limitCorrAv);
       fprintf(fpdnsty,"stepCorr %d\n", stepCorr);
       fprintf(fpdnsty,"nValCorr %d\n", nValCorr);
       fprintf(fpdnsty,"deltaT %lf\n", deltaT);
       real tVal;
       for (n = 1; n <= nValCorr; n++){
     tVal = (n-1)*stepCorr*deltaT;
     fprintf (fpdnsty, "%e\t", tVal);
     int nn = nFunCorr*(n-1);
     for (j = 1; j <= nFunCorr; j ++)
           fprintf (fpdnsty, "%e\t", spacetimeCorrAv[nn + j]);
         fprintf (fpdnsty, "\n");
       }    
       
       countCorrAv = 0.;
       for (j = 1; j <= nFunCorr*nValCorr; j++)
     spacetimeCorrAv[j] = 0.;
     }
   }
  }
}

References cfOrg, cfVal, countCorrAv, deltaT, fpdnsty, indexCorr, limitCorrAv, nAtom, nBuffCorr, NDIM, nFunCorr, nValCorr, region, rx, spacetimeCorr, spacetimeCorrAv, and stepCorr.

EvalUnwrap()

void EvalUnwrap

(

)

Definition at line 27 of file EvalUnwrap.c.

                  {
 int n;
 for (n = 1; n <= nAtom; n++) {
  rxUnwrap[n] = rx[n] + ImageX[n] * region[1];
  ryUnwrap[n] = ry[n] + ImageY[n] * region[2];
  }
}

References ImageX, ImageY, nAtom, region, rx, rxUnwrap, ry, and ryUnwrap.

Referenced by main().

Here is the caller graph for this function:

EvalVrms()

void EvalVrms

(

)

Definition at line 27 of file EvalVrms.c.

               {
  int n;
  VSqr = 0.0;
  VMeanSqr = 0.0;
  VRootMeanSqr = 0.0;
    
  for(n = 1 ; n <= nAtom ; n ++){
  VSqr += Sqr(vx[n]) + Sqr(vy[n]); 
  }
  VMeanSqr = VSqr/nAtom;
  VRootMeanSqr = sqrt(VMeanSqr);
  }

References nAtom, Sqr, VMeanSqr, VRootMeanSqr, VSqr, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

HaltConditionCheck()

bool HaltConditionCheck

(

double

value

)

Definition at line 27 of file Halt.c.

                                      {
  
  if(value <= HaltCondition && value != 0) {
  printf("Halt condition met at step = %d with Vrms = %.16f\n", stepCount, value);
  fprintf(fpresult, "Halt condition met at step = %d with Vrms = %.16f\n", stepCount, value);
  fprintf(fpresult, "Final thermodynamic values:\n");
  fprintf(fpresult, "%0.4lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\n", 
  timeNow, vSum, potEnergy, kinEnergy, totEnergy, uSumPairPerAtom, BondEnergyPerAtom, pressure, virSum);
  return true;       // Signal that the halt condition is met
  }
 return false; // Halt condition not met
}

References BondEnergyPerAtom, fpresult, HaltCondition, kinEnergy, potEnergy, pressure, stepCount, timeNow, totEnergy, uSumPairPerAtom, virSum, and vSum.

Referenced by main().

Here is the caller graph for this function:

Init()

void Init

(

)

Definition at line 31 of file Init.c.

            {
 char dummy[128];
 
 // Always read input parameters
 FILE *fp = fopen("input-data", "r");
 if(!fp) {
  perror("input-data");
  exit(EXIT_FAILURE);
 }
  
  fscanf(fp, "%s %s", mode, inputConfig);  // config type + filename
  fscanf(fp, "%s %s", dummy, solver);
  fscanf(fp, "%s %s %s", dummy, xBoundary, yBoundary);
  fscanf(fp, "%s %d",  dummy, &DampFlag);
  fscanf(fp, "%s %d",  dummy, &freezeAtomType);
  fscanf(fp, "%s %lf", dummy, &rCut);
  fscanf(fp, "%s %lf", dummy, &Kn);
  fscanf(fp, "%s %lf", dummy, &mass);
  fscanf(fp, "%s %lf", dummy, &gamman); 
  fscanf(fp, "%s %lf", dummy, &kappa);
  fscanf(fp, "%s %lf", dummy, &deltaT);
  fscanf(fp, "%s %lf", dummy, &strain);
  fscanf(fp, "%s %lf", dummy, &FyBylx);
  fscanf(fp, "%s %lf", dummy, &fxByfy);
  fscanf(fp, "%s %lf", dummy, &DeltaY);
  fscanf(fp, "%s %lf", dummy, &DeltaX);
  fscanf(fp, "%s %lf", dummy, &HaltCondition);
  fscanf(fp, "%s %d",  dummy, &stepAvg);
  fscanf(fp, "%s %d",  dummy, &stepEquil);
  fscanf(fp, "%s %d",  dummy, &stepLimit);
  fscanf(fp, "%s %d",  dummy, &stepDump);
  fscanf(fp, "%s %d",  dummy, &stepTraj);
  fscanf(fp, "%s %d",  dummy, &limitCorrAv);
  fscanf(fp, "%s %d",  dummy, &nBuffCorr);
  fscanf(fp, "%s %d",  dummy, &nFunCorr);
  fscanf(fp, "%s %d",  dummy, &nValCorr);
  fscanf(fp, "%s %d",  dummy, &stepCorr);
  fscanf(fp, "%s %d",  dummy, &limitAcfAv);
  fscanf(fp, "%s %d",  dummy, &nBuffAcf);
  fscanf(fp, "%s %d",  dummy, &nValAcf);
  fscanf(fp, "%s %d",  dummy, &stepAcf);
  fscanf(fp, "%s %lf", dummy, &rangeRdf);
  fscanf(fp, "%s %d",  dummy, &limitRdf);
  fscanf(fp, "%s %d",  dummy, &sizeHistRdf);
  fscanf(fp, "%s %d",  dummy, &stepRdf);
  fclose(fp);
  
  int useBinaryRestart = 0;
  if(strcmp(mode, "read_restart") == 0) {
   useBinaryRestart = 1;
  } else if (strcmp(mode, "read_data") != 0) {
   fprintf(stderr, "ERROR: First line of input-data must be 'read_data' or 'read_restart'\n");
   exit(0);
  }
    
  //Conditionally read binary config
  if(useBinaryRestart) {
   ReadBinaryRestart(inputConfig);  // uses global prefix + config file
   printf(">>> Binary restart loaded from %s <<<\n", inputConfig);
   printf(">>> Restarting simulation from time = %.8lf <<<\n", timeNow);
   return; //Exiting from Init() from here 
  }
  
  FILE *fpSTATE;
  if((fpSTATE = fopen(inputConfig,"r"))==NULL){
  printf("Error occurred: Could not open STATE file\n Exiting now..\n");
  exit(0);
  }
 
  if(fscanf(fpSTATE, "%s %lf", dummy, &timeNow) != 2 || strcmp(dummy, "timeNow") != 0) {
   fprintf(stderr, "ERROR [%s:%d:%s]: Expected 'timeNow <value>' as the first line in the config file.\n",
   __FILE__, __LINE__, __func__);
   exit(EXIT_FAILURE);
  }
 
  if(timeNow == 0.0) {
   printf(">>> Running from time = 0.0: Beginning of the simulation\n");
   stepCount = 0;
  } 
  
  fscanf(fpSTATE, "%s %d",  dummy, &nAtom);
  fscanf(fpSTATE, "%s %d",  dummy, &nBond);
  fscanf(fpSTATE, "%s %d",  dummy, &nAtomType);
  fscanf(fpSTATE, "%s %d",  dummy, &nBondType);
  fscanf(fpSTATE, "%s %lf", dummy, &region[1]);
  fscanf(fpSTATE, "%s %lf", dummy, &region[2]);
  
  if(timeNow == 0.0)  region[2] *= 1.5; //Remove this when put on GitHub
 
  density = nAtom/(region[1]*region[2]); 
  cells[1] = region[1] / rCut;
  cells[2] = region[2] / rCut;
  cellList = (int*)malloc((nAtom + cells[1] * cells[2] + 1) * sizeof(int));
  regionH[1] = 0.5*region[1];
  regionH[2] = 0.5*region[2];
 
  //strain information
  strainRate = strain/deltaT;
  shearDisplacement = strain * region[2];
  shearVelocity = strainRate * region[2];
  int n;
 
  rx = (double*)malloc((nAtom + 1) * sizeof(double));
  ry = (double*)malloc((nAtom + 1) * sizeof(double));
  vx = (double*)malloc((nAtom + 1) * sizeof(double));
  vy = (double*)malloc((nAtom + 1) * sizeof(double));
  ax = (double*)malloc((nAtom + 1) * sizeof(double));
  ay = (double*)malloc((nAtom + 1) * sizeof(double));
  fx = (double*)malloc((nAtom + 1) * sizeof(double));
  fy = (double*)malloc((nAtom + 1) * sizeof(double));
  fax = (double*)malloc((nAtom + 1) * sizeof(double));
  fay = (double*)malloc((nAtom + 1) * sizeof(double));
  atomID = (int*)malloc((nAtom+1) * sizeof(int));
  atomType = (int*)malloc((nAtom+1) * sizeof(int));
  atomRadius = (double*)malloc((nAtom + 1) * sizeof(double));
  atomMass = (double*)malloc((nAtom + 1) * sizeof(double));
  speed = (double*)malloc((nAtom + 1) * sizeof(double));
  discDragx = (double*)malloc((nAtom + 1) * sizeof(double)); 
  discDragy = (double*)malloc((nAtom + 1) * sizeof(double)); 
  molID = (int*)malloc((nAtom+1) * sizeof(int));
    
  BondID = (int*)malloc((nBond+1)*sizeof(int));
  BondType = (int*)malloc((nBond+1)*sizeof(int));
  atom1 = (int*)malloc((nBond+1)*sizeof(int));
  atom2 = (int*)malloc((nBond+1)*sizeof(int)); 
  kb = (double*)malloc((nBond+1)*sizeof(double));
  ro = (double*)malloc((nBond+1)*sizeof(double));
  BondEnergy = (double*)malloc((nBond+1)*sizeof(double));
  BondLength =(double*)malloc((nBond+1)*sizeof(double));
  nodeDragx = (double*)malloc((nAtom + 1) * sizeof(double));
  nodeDragy = (double*)malloc((nAtom + 1) * sizeof(double));
  rxUnwrap = (double*)malloc((nAtom + 1) * sizeof(double));
  ryUnwrap = (double*)malloc((nAtom + 1) * sizeof(double));
  ImageX = (int*)malloc((nAtom+1) * sizeof(int));
  ImageY = (int*)malloc((nAtom+1) * sizeof(int));
 
 
  for(n = 1; n <= nAtom; n ++){
   atomMass[n] = mass;
  }
 
  fscanf(fpSTATE, "%s\n", dummy);
  for(n = 1; n <= nAtom; n ++)
   fscanf(fpSTATE, "%d %d %d %lf %lf %lf %lf %lf\n", &atomID[n], &molID[n], &atomType[n], &atomRadius[n], &rx[n], &ry[n], &vx[n], &vy[n]);
   
 
  fscanf(fpSTATE, "%s\n", dummy); 
  for(n=1; n<=nBond; n++)
   fscanf(fpSTATE, "%d %d %d %d %lf %lf\n", &BondID[n], &BondType[n], &atom1[n], &atom2[n], &kb[n], &ro[n]);
 
  fclose(fpSTATE);
 
 //2D-List of bonded atoms. This is used to remove pair interaction
 //calculation for the bonded atoms
   isBonded = (int**)malloc((nAtom + 1) * sizeof(int*));
  for (int i = 0; i <= nAtom; i++) {
    isBonded[i] = (int*)malloc((nAtom + 1) * sizeof(int));
    for (int j = 0; j <= nAtom; j++) {
        isBonded[i][j] = 0;
    }
  }
  
  for (n = 1; n <= nBond; n++) {
    int i = atom1[n];
    int j = atom2[n];
    isBonded[i][j] = 1;
    isBonded[j][i] = 1; // symmetric
}
 
 // List the interface atoms
 nAtomInterface = 0;
 nAtomBlock = 0;
 nDiscInterface = 0;
 bigDiameter = 2.8;
 InterfaceWidth = 5.0 * bigDiameter;
 
 for(n = 1; n <= nAtom; n++){
  if(fabs(ry[n]) < InterfaceWidth){
  nAtomInterface++;
  }
  if(molID[n] == 2){
  nAtomBlock++;
  } 
  if(atomRadius[n] != 0.0){
  nDiscInterface++;
  } } 
 
  
  int BondPairInteract;
  BondPairInteract = 0;
  int m;
  if(BondPairInteract == 1){
   atomIDInterface =  (int *)malloc((nAtomInterface+1)*sizeof(int));
   m = 1;
   for(n=1; n<=nAtom; n++){
    if(fabs(ry[n]) < InterfaceWidth){
    atomIDInterface[m] = atomID[n]; 
    m++; 
  } } }
  else if(BondPairInteract == 0){
   nAtomInterface = nDiscInterface;
   atomIDInterface =  (int *)malloc((nAtomInterface+1)*sizeof(int));
   m = 1;
   for(n=1; n<=nAtom; n++){
    if(atomRadius[n] != 0.0){
    atomIDInterface[m] = atomID[n]; 
    m++; 
  } } }
 
  nPairTotal = 0.5 * nAtomInterface * (nAtomInterface-1);
  PairID = (int*)malloc((nPairTotal+1) * sizeof(int));
  Pairatom1 = (int*)malloc((nPairTotal+1) * sizeof(int));
  Pairatom2 = (int*)malloc((nPairTotal+1) * sizeof(int));
  PairXij = (double*)malloc((nPairTotal+1) * sizeof(double));
  PairYij = (double*)malloc((nPairTotal+1) * sizeof(double));
  
  fprintf(fpresult, "------------------------------------\n");
  fprintf(fpresult, "-------PARAMETERS-----------\n");
  fprintf(fpresult, "------------------------------------\n");
  fprintf(fpresult, "nAtom\t\t\t%d\n",  nAtom);
  fprintf(fpresult, "nBond\t\t\t%d\n",  nBond);
  fprintf(fpresult, "nAtomType\t\t%d\n", nAtomType);
  fprintf(fpresult, "nBondType\t\t%d\n",  nBondType);
  fprintf(fpresult, "nAtomBlock\t\t%d\n",  nAtomBlock);
  fprintf(fpresult, "nAtomInterface\t\t%d\n",  nAtomInterface);
  fprintf(fpresult, "nDiscInterface\t\t%d\n",  nDiscInterface);
  fprintf(fpresult, "mass\t\t\t%0.6g\n", mass);
  fprintf(fpresult, "gamman\t\t\t%0.6g\n", gamman);
  fprintf(fpresult, "strain\t\t\t%0.6g\n", strain);
  fprintf(fpresult, "strainRate\t\t%0.6g\n", strainRate);
  fprintf(fpresult, "FyBylx\t\t\t%0.6g\n", FyBylx);
  fprintf(fpresult, "fxByfy\t\t\t%0.6g\n", fxByfy);
  fprintf(fpresult, "DeltaY\t\t\t%0.6g\n", DeltaY);
  fprintf(fpresult, "DeltaX\t\t\t%0.6g\n", DeltaX);
  fprintf(fpresult, "HaltCondition\t\t%0.6g\n", HaltCondition);
  fprintf(fpresult, "kappa\t\t\t%g\n", kappa);
  fprintf(fpresult, "density\t\t\t%g\n", density);
  fprintf(fpresult, "rCut\t\t\t%g\n", rCut);
  fprintf(fpresult, "deltaT\t\t\t%g\n", deltaT);
  fprintf(fpresult, "stepEquil\t\t%d\n",  stepEquil);
  fprintf(fpresult, "stepLimit\t\t%d\n",  stepLimit);
  fprintf(fpresult, "region[1]\t\t%0.16lf\n", region[1]);
  fprintf(fpresult, "region[2]\t\t%0.16lf\n", region[2]);
  fprintf(fpresult, "cells[1]\t\t%d\n",  cells[1]);
  fprintf(fpresult, "cells[2]\t\t%d\n",  cells[2]);
  fprintf(fpresult, "solver\t\t\t%s\n",  solver);
  fprintf(fpresult, "boundary\t\t%s %s\n", xBoundary, yBoundary);
  fprintf(fpresult, "DampFlag\t\t%d\n",  DampFlag);
  fprintf(fpresult, "------------------------------------\n");
  fprintf(fpresult, "#TimeNow TotalMomentum PotEngyPerAtom KinEngyPerAtom TotEngyPerAtom PairEnergyPerAtom BondEnergyPerAtom  Press VirialSum\n");
  fprintf(fpvrms, "#timeNow\tVrms \n");
  fprintf(fpcom, "#timeNow\tComX\tComY\n");
  fprintf(fpforce, "#timeNow\tforceSumxAtomType1\tforceSumyAtomType1\tforceSumxAtomType2\tforceSumyAtomType2\tforceSumxAtomType3\tforceSumyAtomType3\tforceSumxAtomType4\tforceSumyAtomType4\tforceSumxAtomType5\tforceSumyAtomType5\tforceSumxExtern\tforceSumyExtern\n");
 
/* //Uncomment the following as per your acquirement
    fprintf(fpstress, "strain           %lf\n", strain);
    fprintf(fpstress, "region[1]        %lf\n", region[1]);
    fprintf(fpstress, "region[2]        %lf\n", region[2]);
    fprintf(fpstress, "#timeNow virSumxx virSumyy virSumxy pressure\n");
    fprintf(fpmomentum, "#timeNow Px Py\n");
*/
  
   if((strcmp(xBoundary, "p") != 0 && strcmp(xBoundary, "r") != 0) ||
    (strcmp(yBoundary, "p") != 0 && strcmp(yBoundary, "r") != 0)) {
    fprintf(fpresult, "Error: Invalid boundary value detected: '%s %s'. Only 'p' or 'r' are allowed.\n", xBoundary, yBoundary);
    exit(EXIT_FAILURE); // Exit with failure status
  }
 
}

References atom1, atom2, atomID, atomIDInterface, atomMass, atomRadius, atomType, ax, ay, bigDiameter, BondEnergy, BondID, BondLength, BondType, cellList, cells, DampFlag, deltaT, DeltaX, DeltaY, density, discDragx, discDragy, fax, fay, fpcom, fpforce, fpresult, fpvrms, freezeAtomType, fx, fxByfy, fy, FyBylx, gamman, HaltCondition, ImageX, ImageY, inputConfig, InterfaceWidth, isBonded, kappa, kb, Kn, limitAcfAv, limitCorrAv, limitRdf, mass, mode, molID, nAtom, nAtomBlock, nAtomInterface, nAtomType, nBond, nBondType, nBuffAcf, nBuffCorr, nDiscInterface, nFunCorr, nodeDragx, nodeDragy, nPairTotal, nValAcf, nValCorr, Pairatom1, Pairatom2, PairID, PairXij, PairYij, rangeRdf, rCut, ReadBinaryRestart(), region, regionH, ro, rx, rxUnwrap, ry, ryUnwrap, shearDisplacement, shearVelocity, sizeHistRdf, solver, speed, stepAcf, stepAvg, stepCorr, stepCount, stepDump, stepEquil, stepLimit, stepRdf, stepTraj, strain, strainRate, timeNow, vx, vy, xBoundary, and yBoundary.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

main()

int main	(	int	argc,
		char **	argv )

Definition at line 43 of file main.c.

                                {
 time_t t1 = 0, t2;
 if (argc < 2) {
 fprintf(stderr, "Usage: %s <output_prefix>\n", argv[0]);
 return 1;
 }
 
 int prefix_size = snprintf(NULL, 0, "../output/%s", argv[1]) + 1; // +1 for the null terminator
 prefix = malloc(prefix_size);
 if(prefix == NULL) {
  fprintf(stderr, "Memory allocation failed\n");
 return 1;
 }
 
   // Write the formatted string into the allocated space
  snprintf(prefix, prefix_size, "../output/%s", argv[1]);
  sprintf(result, "%s.result", prefix);
  fpresult = fopen(result, "w");
  sprintf(xyz, "%s.xyz", prefix);
  fpxyz = fopen(xyz, "w");
  sprintf(vrms, "%s.vrms", prefix);
  fpvrms = fopen(vrms, "w");
  sprintf(bond, "%s.bond", prefix);
  fpbond = fopen(bond, "w");
  sprintf(com, "%s.com", prefix);
  fpcom = fopen(com, "w");
  sprintf(pair, "%s.pair", prefix);
  fppair = fopen(pair, "w");
  sprintf(force, "%s.force", prefix);
  fpforce = fopen(force, "w");
  
  /* //Uncomment the following as per your acquirement 
  sprintf(dnsty, "%s.curr-dnsty", prefix);
  fpdnsty = fopen(dnsty, "w");
  sprintf(visc, "%s.viscosity", prefix);
  fpvisc = fopen(visc, "w");
  sprintf(rdf, "%s.rdf", prefix);
  fprdf = fopen(rdf, "w");
  sprintf(stress, "%s.stress", prefix);
  fpstress = fopen(stress, "w");
  sprintf(momentum, "%s.momentum", prefix);
  fpmomentum = fopen(momentum, "w");
  */  
 
  Init();
  SetupJob();
  t1 = time(NULL);
  moreCycles = 1;
  if(stepCount >= 0) {
   if (timeNow == 0.0) {
    printf(">>> Run type: Fresh simulation <<<\n");
    DisplaceAtoms();
    ComputePairForce(1);
    ComputeBondForce();
    ApplyForce();
    } else {
    printf(">>> Run type: Restart simulation <<<\n");
   }
   DumpBonds();
   DumpPairs();
   Trajectory();
   EvalUnwrap();
   ApplyBoundaryCond();
   EvalProps();
   EvalVrms();
   EvalCom();
   PrintVrms();
   PrintCom();
   PrintSummary();
   PrintForceSum();
   }
 
//Here starts the main loop of the program 
  while(moreCycles){
   if(stepLimit == 0){
   printf("Error occured: stepLimit must be > 0\n");
   printf("Exiting now ...\n");
   exit(0);
   }
 
   stepCount ++;
   timeNow += deltaT ; //stepCount * deltaT; //for adaptive step size: timeNow += deltaT
 
   VelocityVerletStep(1);
   EvalUnwrap();
   ApplyBoundaryCond();
   ComputePairForce(1);
   ComputeBondForce();
   ApplyForce();
   VelocityVerletStep(2);
   ApplyBoundaryCond();
   EvalProps();
   EvalVrms();
   EvalCom();
   if(stepCount % stepAvg == 0){
    PrintSummary();
    PrintVrms();
    PrintCom();
    PrintForceSum();
    }
   if(stepCount % stepTraj == 0){
    Trajectory();
    DumpBonds();
    DumpPairs();
    }
   if(stepCount % stepDump == 0){
    DumpRestart();     // Save the current state for input
    DumpState();       // Save the current state for config
    WriteBinaryRestart();
   }
   if(HaltConditionCheck(VRootMeanSqr)) {
    DumpRestart();     // Save the current state for input
    DumpState();       // Save the current state for config
    WriteBinaryRestart();
    break;  // Exit the loop when the halt condition is met
    }
    
    moreCycles ++;
    if(moreCycles >= stepLimit)
    moreCycles = 0;
  }
 
 
  t2 = time(NULL);
  fprintf(fpresult, "#Execution time %lf secs\n", difftime(t2,t1));
  fprintf(fpresult, "#Execution speed %lf steps per secs\n", stepLimit/difftime(t2,t1));
  printf(">>> Simulation run completed <<<\n");
  printf(">>> Execution time %lf secs <<<\n", difftime(t2,t1));
  printf(">>> Execution speed %lf steps per secs <<< \n", stepLimit/difftime(t2,t1));
 
  fclose(fpresult);
  fclose(fpxyz);
  fclose(fpvrms);
  fclose(fpbond);
  fclose(fppair);
  fclose(fpcom);
  fclose(fpforce);
 
/*//Uncomment the following as per your acquirement    
  fclose(fpdnsty);
  fclose(fpvisc);
  fclose(fprdf);
  fclose(fpstress);
  fclose(fpmomentum);
*/
 
  free(prefix);
  Close();
  return 0;
} 

References ApplyBoundaryCond(), ApplyForce(), bond, Close(), com, ComputeBondForce(), ComputePairForce(), deltaT, DisplaceAtoms(), DumpBonds(), DumpPairs(), DumpRestart(), DumpState(), EvalCom(), EvalProps(), EvalUnwrap(), EvalVrms(), force, fpbond, fpcom, fpforce, fppair, fpresult, fpvrms, fpxyz, HaltConditionCheck(), Init(), moreCycles, pair, prefix, PrintCom(), PrintForceSum(), PrintSummary(), PrintVrms(), result, SetupJob(), stepAvg, stepCount, stepDump, stepLimit, stepTraj, timeNow, Trajectory(), VelocityVerletStep(), vrms, VRootMeanSqr, WriteBinaryRestart(), and xyz.

Here is the call graph for this function:

PrintCom()

void PrintCom

(

)

Definition at line 28 of file PrintCom.c.

               {
 fprintf(fpcom, "%0.4lf\t%0.16lf\t%0.16lf\n", timeNow, ComX, ComY);
 fflush(fpcom);
 }

References ComX, ComY, fpcom, and timeNow.

Referenced by main().

Here is the caller graph for this function:

PrintForceSum()

void PrintForceSum

(

)

Definition at line 28 of file PrintForceSum.c.

                    {
 int n;
 double forceSumxAtomType1, forceSumxAtomType2, forceSumxAtomType3, forceSumxAtomType4, forceSumxAtomType5;
 double forceSumyAtomType1, forceSumyAtomType2, forceSumyAtomType3, forceSumyAtomType4, forceSumyAtomType5;
 
 forceSumxAtomType1 = 0.0;  forceSumyAtomType1 = 0.0; 
 forceSumxAtomType2 = 0.0;  forceSumyAtomType2 = 0.0; 
 forceSumxAtomType3 = 0.0;  forceSumyAtomType3 = 0.0; 
 forceSumxAtomType4 = 0.0;  forceSumyAtomType4 = 0.0; 
 forceSumxAtomType5 = 0.0;  forceSumyAtomType5 = 0.0; 
 
 
 for(n = 1; n <= nAtom; n++){
 if(atomType[n] == 1){
  forceSumxAtomType1 += fx[n];
  forceSumyAtomType1  += fy[n];
  } else if(atomType[n] == 2){
  forceSumxAtomType2 += fx[n];
  forceSumyAtomType2  += fy[n];
  } else if(atomType[n] == 3){
  forceSumxAtomType3 += fx[n];
  forceSumyAtomType3  += fy[n];
  } else if(atomType[n] == 4){
  forceSumxAtomType4 += fx[n];
  forceSumyAtomType4 += fy[n];
  } else if(atomType[n] == 5){
  forceSumxAtomType5 += fx[n];
  forceSumyAtomType5 += fy[n];
  }
 } 
 
 fprintf(fpforce, "%0.4lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16f\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16f\t%0.16f\n", timeNow, 
 forceSumxAtomType1, forceSumyAtomType1, 
 forceSumxAtomType2, forceSumyAtomType2, 
 forceSumxAtomType3, forceSumyAtomType3, 
 forceSumxAtomType4, forceSumyAtomType4, 
 forceSumxAtomType5, forceSumyAtomType5, 
 forceSumxExtern, forceSumyExtern);
 fflush(fpforce);
 }

References atomType, forceSumxExtern, forceSumyExtern, fpforce, fx, fy, nAtom, and timeNow.

Referenced by main().

Here is the caller graph for this function:

PrintMomentum()

void PrintMomentum

(

)

Definition at line 25 of file PrintMomentum.c.

                    {
 fprintf(fpmomentum, "%0.4lf\t%0.16lf\t%0.16lf\n", timeNow, vSumX, vSumY);
 fflush(fpmomentum);
}

References fpmomentum, timeNow, vSumX, and vSumY.

PrintStress()

void PrintStress

(

)

Definition at line 25 of file PrintStress.c.

                  {
  fprintf(fpstress, "%0.4lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\n", timeNow, virSumxx, virSumyy, virSumxy, pressure);
  fflush(fpstress);
}

References fpstress, pressure, timeNow, virSumxx, virSumxy, and virSumyy.

PrintSummary()

void PrintSummary

(

)

Definition at line 25 of file PrintSummary.c.

                   {
 fprintf(fpresult, "%0.4lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\t%0.16lf\n", 
  timeNow, vSum, potEnergy, kinEnergy, totEnergy, uSumPairPerAtom, BondEnergyPerAtom, pressure, virSum);
  fflush(fpresult);
}

References BondEnergyPerAtom, fpresult, kinEnergy, potEnergy, pressure, timeNow, totEnergy, uSumPairPerAtom, virSum, and vSum.

Referenced by main().

Here is the caller graph for this function:

PrintVrms()

void PrintVrms

(

)

Definition at line 27 of file PrintVrms.c.

                {
 fprintf(fpvrms, "%0.4lf\t%0.16lf\n", timeNow, VRootMeanSqr);
 fflush(fpvrms);
}

References fpvrms, timeNow, and VRootMeanSqr.

Referenced by main().

Here is the caller graph for this function:

SetupJob()

void SetupJob

(

)

Definition at line 27 of file SetupJob.c.

               {
  AllocArrays();
  AccumProps(0);
  InitVacf();
  // INITIALISE SPACETIME CORRELATIONS
  int n;
  for (n = 1; n <= nBuffCorr; n++)
    indexCorr[n] = -(n - 1)*nValCorr/nBuffCorr;
 
  countCorrAv = 0.;
 
  for (n = 1; n <= nFunCorr*nValCorr; n++)
    spacetimeCorrAv[n] = 0.;
 
  //RDF
  countRdf = 0;
}

References AccumProps(), AllocArrays(), countCorrAv, countRdf, indexCorr, InitVacf(), nBuffCorr, nFunCorr, nValCorr, and spacetimeCorrAv.

Referenced by main().

Here is the call graph for this function:

Here is the caller graph for this function:

Trajectory()

void Trajectory

(

)

Definition at line 25 of file Trajectory.c.

                 {
 int n;
 //Trajectory file in LAMMPS dump format for OVITO visualization
 fprintf(fpxyz, "ITEM: TIMESTEP\n");
 fprintf(fpxyz, "%lf\n",timeNow);
 fprintf(fpxyz, "ITEM: NUMBER OF ATOMS\n");
 fprintf(fpxyz, "%d\n",nAtom);
 fprintf(fpxyz, "ITEM: BOX BOUNDS pp ff pp\n");
 fprintf(fpxyz, "%lf %lf xlo xhi\n", -regionH[1], regionH[1]);
 fprintf(fpxyz, "%lf %lf ylo yhi\n", -regionH[2], regionH[2]);
 fprintf(fpxyz, "%lf %lf zlo zhi\n", -0.1, 0.1);
 fprintf(fpxyz, "ITEM: ATOMS id mol type radius x y vx vy fx fy\n");
 for(n=1; n<=nAtom; n++)
  fprintf(fpxyz, "%d\t %d\t %d\t %0.2lf\t %0.16lf\t %0.16lf\t %0.16lf\t %0.16lf\t %0.16lf\t %0.16lf\n", 
  atomID[n], molID[n], atomType[n], atomRadius[n], rx[n], ry[n], vx[n], vy[n], fx[n], fy[n]);
}

References atomID, atomRadius, atomType, fpxyz, fx, fy, molID, nAtom, regionH, rx, ry, timeNow, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

VelocityVerletStep()

void VelocityVerletStep

(

int

icode

)

Definition at line 26 of file VelocityVerletStep.c.

                                  {
int n;
double atomMassi;
   
 if(icode == 1){
 for (n= 1; n <= nAtom; n++) {
  if(atomType[n] != freezeAtomType){
  atomMassi = 1./atomMass[n];
  ax[n] = fx[n] * atomMassi;    ay[n] = fy[n] * atomMassi;    
  vx[n] += ax[n] * 0.5 * deltaT; 
  vy[n] += ay[n] * 0.5 * deltaT;
  rx[n] += vx[n] * deltaT; 
  ry[n] += vy[n] * deltaT;
  } 
  //Calculating the image flags here
  if (rx[n] >= regionH[1]) {
   rx[n] -= region[1];
   ImageX[n]++;
   } else if (rx[n] < -regionH[1]) {
   rx[n] += region[1];
  ImageX[n]--;
  }
  if (ry[n] >= regionH[2]) {
   ry[n] -= region[2];
   ImageY[n]++;
   } else if (ry[n] < -regionH[2]) {
   ry[n] += region[2];
   ImageY[n]--;
   } } }
  else if(icode == 2){
  for(n = 1; n <= nAtom; n++) {
  if(atomType[n] != freezeAtomType){
   atomMassi = 1./atomMass[n];
   ax[n] = fx[n] * atomMassi;    ay[n] = fy[n] * atomMassi;    
   vx[n] += ax[n] * 0.5 * deltaT; 
   vy[n] += ay[n] * 0.5 * deltaT;
} } } }

References atomMass, atomType, ax, ay, deltaT, freezeAtomType, fx, fy, ImageX, ImageY, nAtom, region, regionH, rx, ry, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

WriteBinaryRestart()

void WriteBinaryRestart

(

)

Definition at line 60 of file WriteRestartBinary.c.

                          {
 RestartHeader hdr = {
 .magic = "LAMINA",
 .version = 1.0,
 .timeNow = timeNow,
 .nAtom = nAtom,
 .nBond = nBond,
 .nAtomType = nAtomType,
 .nBondType = nBondType,
 .regionX = region[1],
 .regionY = region[2],
 .nAtomInterface = nAtomInterface,
 .nAtomBlock = nAtomBlock,
 .nDiscInterface = nDiscInterface,
 .bigDiameter = bigDiameter,
 .InterfaceWidth = InterfaceWidth,
 .nPairActive = nPairActive,
 .nPairTotal = nPairTotal,
 .uSumPair = uSumPair,
 .virSumPair = virSumPair,
 .virSumPairxx = virSumPairxx,
 .virSumPairyy = virSumPairyy,
 .virSumPairxy = virSumPairxy,
 .TotalBondEnergy = TotalBondEnergy,
 .virSumBond = virSumBond,
 .virSumBondxx = virSumBondxx,
 .virSumBondyy = virSumBondyy,
 .virSumBondxy = virSumBondxy,
 .stepCount = stepCount,
 .forceSumxExtern = forceSumxExtern, 
 .forceSumyExtern = forceSumyExtern
 };
 
 char DUMP[256];
 FILE *fp;
 sprintf(DUMP, "%s.bin", prefix);  // Produces e.g. "../output/test.bin"
 fp = fopen(DUMP, "wb");
 if (!fp) {
 fprintf(stderr, "Error opening binary restart file %s for writing\n", DUMP);
 exit(EXIT_FAILURE);
 }
 
 //Here we are writing the data to binary file
 fwrite(&hdr, sizeof(RestartHeader), 1, fp);
 fwrite(&atomID[1], sizeof(int), nAtom, fp); 
 fwrite(&molID[1], sizeof(int), nAtom, fp); 
 fwrite(&atomType[1], sizeof(int), nAtom, fp); 
 fwrite(&atomRadius[1], sizeof(double), nAtom, fp); 
 fwrite(&rx[1], sizeof(double), nAtom, fp); 
 fwrite(&ry[1], sizeof(double), nAtom, fp); 
 fwrite(&vx[1], sizeof(double), nAtom, fp); 
 fwrite(&vy[1], sizeof(double), nAtom, fp); 
 fwrite(&ax[1], sizeof(double), nAtom, fp); 
 fwrite(&ay[1], sizeof(double), nAtom, fp); 
 fwrite(&fx[1], sizeof(double), nAtom, fp); 
 fwrite(&fy[1], sizeof(double), nAtom, fp); 
 fwrite(&atomMass[1], sizeof(double), nAtom, fp);
 fwrite(&discDragx[1], sizeof(double), nAtom, fp);
 fwrite(&discDragy[1], sizeof(double), nAtom, fp);
 fwrite(&atomIDInterface[1], sizeof(int), nAtomInterface, fp); 
 
 fwrite(&BondID[1], sizeof(int), nBond, fp);
 fwrite(&BondType[1], sizeof(int), nBond, fp);
 fwrite(&atom1[1], sizeof(int), nBond, fp);
 fwrite(&atom2[1], sizeof(int), nBond, fp);
 fwrite(&kb[1], sizeof(double), nBond, fp);
 fwrite(&ro[1], sizeof(double), nBond, fp); 
 fwrite(&BondEnergy[1], sizeof(double), nBond, fp);
 fwrite(&BondLength[1], sizeof(double), nBond, fp);    
 fwrite(&nodeDragx[1], sizeof(double), nAtom, fp);
 fwrite(&nodeDragy[1], sizeof(double), nAtom, fp);
 fwrite(&rxUnwrap[1], sizeof(double), nAtom, fp);
 fwrite(&ryUnwrap[1], sizeof(double), nAtom, fp);
 fwrite(&ImageX[1], sizeof(int), nAtom, fp); 
 fwrite(&ImageY[1], sizeof(int), nAtom, fp); 
 
 fwrite(&PairID[1], sizeof(int), nPairActive, fp);
 fwrite(&Pairatom1[1], sizeof(int), nPairActive, fp);
 fwrite(&Pairatom2[1], sizeof(int), nPairActive, fp);
 fwrite(&PairXij[1], sizeof(double), nPairActive, fp);
 fwrite(&PairYij[1], sizeof(double), nPairActive, fp);
 
 fclose(fp);
}

References atom1, atom2, atomID, atomIDInterface, atomMass, atomRadius, atomType, ax, ay, bigDiameter, BondEnergy, BondID, BondLength, BondType, discDragx, discDragy, forceSumxExtern, forceSumyExtern, fx, fy, ImageX, ImageY, InterfaceWidth, kb, molID, nAtom, nAtomBlock, nAtomInterface, nAtomType, nBond, nBondType, nDiscInterface, nodeDragx, nodeDragy, nPairActive, nPairTotal, Pairatom1, Pairatom2, PairID, PairXij, PairYij, prefix, region, ro, rx, rxUnwrap, ry, ryUnwrap, stepCount, timeNow, TotalBondEnergy, uSumPair, virSumBond, virSumBondxx, virSumBondxy, virSumBondyy, virSumPair, virSumPairxx, virSumPairxy, virSumPairyy, vx, and vy.

Referenced by main().

Here is the caller graph for this function:

Variable Documentation

prefix

char* prefix = NULL

Definition at line 13 of file main.c.

Referenced by DumpRestart(), DumpState(), main(), and WriteBinaryRestart().