Lemina
A molecular dynamics package for network, granular material and point particles with a range of interaction potential.
 
Loading...
Searching...
No Matches
Init.c
Go to the documentation of this file.
1/*
2 * This file is part of Lamina.
3 *
4 * Lamina is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * Lamina is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with Lamina. If not, see <https://www.gnu.org/licenses/>.
16
17 Copyright (C) 2025 Harish Charan, University of Durham, UK
18
19 */
20
21
22#include<stdio.h>
23#include<math.h>
24#include<string.h>
25#include<stdlib.h>
26#include <time.h>
27#include"global.h"
28
29void ReadBinaryRestart(const char *filename);
30
31void Init() {
32 char dummy[128];
33
34 // Always read input parameters
35 FILE *fp = fopen("input-data", "r");
36 if(!fp) {
37 perror("input-data");
38 exit(EXIT_FAILURE);
39 }
40
41 fscanf(fp, "%s %s", mode, inputConfig); // config type + filename
42 fscanf(fp, "%s %s", dummy, solver);
43 fscanf(fp, "%s %s %s", dummy, xBoundary, yBoundary);
44 fscanf(fp, "%s %d", dummy, &DampFlag);
45 fscanf(fp, "%s %d", dummy, &freezeAtomType);
46 fscanf(fp, "%s %lf", dummy, &rCut);
47 fscanf(fp, "%s %lf", dummy, &Kn);
48 fscanf(fp, "%s %lf", dummy, &mass);
49 fscanf(fp, "%s %lf", dummy, &gamman);
50 fscanf(fp, "%s %lf", dummy, &kappa);
51 fscanf(fp, "%s %lf", dummy, &deltaT);
52 fscanf(fp, "%s %lf", dummy, &strain);
53 fscanf(fp, "%s %lf", dummy, &FyBylx);
54 fscanf(fp, "%s %lf", dummy, &fxByfy);
55 fscanf(fp, "%s %lf", dummy, &DeltaY);
56 fscanf(fp, "%s %lf", dummy, &DeltaX);
57 fscanf(fp, "%s %lf", dummy, &HaltCondition);
58 fscanf(fp, "%s %d", dummy, &stepAvg);
59 fscanf(fp, "%s %d", dummy, &stepEquil);
60 fscanf(fp, "%s %d", dummy, &stepLimit);
61 fscanf(fp, "%s %d", dummy, &stepDump);
62 fscanf(fp, "%s %d", dummy, &stepTraj);
63 fscanf(fp, "%s %d", dummy, &limitCorrAv);
64 fscanf(fp, "%s %d", dummy, &nBuffCorr);
65 fscanf(fp, "%s %d", dummy, &nFunCorr);
66 fscanf(fp, "%s %d", dummy, &nValCorr);
67 fscanf(fp, "%s %d", dummy, &stepCorr);
68 fscanf(fp, "%s %d", dummy, &limitAcfAv);
69 fscanf(fp, "%s %d", dummy, &nBuffAcf);
70 fscanf(fp, "%s %d", dummy, &nValAcf);
71 fscanf(fp, "%s %d", dummy, &stepAcf);
72 fscanf(fp, "%s %lf", dummy, &rangeRdf);
73 fscanf(fp, "%s %d", dummy, &limitRdf);
74 fscanf(fp, "%s %d", dummy, &sizeHistRdf);
75 fscanf(fp, "%s %d", dummy, &stepRdf);
76 fclose(fp);
77
78 int useBinaryRestart = 0;
79 if(strcmp(mode, "read_restart") == 0) {
80 useBinaryRestart = 1;
81 } else if (strcmp(mode, "read_data") != 0) {
82 fprintf(stderr, "ERROR: First line of input-data must be 'read_data' or 'read_restart'\n");
83 exit(0);
84 }
85
86 //Conditionally read binary config
87 if(useBinaryRestart) {
88 ReadBinaryRestart(inputConfig); // uses global prefix + config file
89 printf(">>> Binary restart loaded from %s <<<\n", inputConfig);
90 printf(">>> Restarting simulation from time = %.8lf <<<\n", timeNow);
91 return; //Exiting from Init() from here
92 }
93
94 FILE *fpSTATE;
95 if((fpSTATE = fopen(inputConfig,"r"))==NULL){
96 printf("Error occurred: Could not open STATE file\n Exiting now..\n");
97 exit(0);
98 }
99
100 if(fscanf(fpSTATE, "%s %lf", dummy, &timeNow) != 2 || strcmp(dummy, "timeNow") != 0) {
101 fprintf(stderr, "ERROR [%s:%d:%s]: Expected 'timeNow <value>' as the first line in the config file.\n",
102 __FILE__, __LINE__, __func__);
103 exit(EXIT_FAILURE);
104 }
105
106 if(timeNow == 0.0) {
107 printf(">>> Running from time = 0.0: Beginning of the simulation\n");
108 stepCount = 0;
109 }
110
111 fscanf(fpSTATE, "%s %d", dummy, &nAtom);
112 fscanf(fpSTATE, "%s %d", dummy, &nBond);
113 fscanf(fpSTATE, "%s %d", dummy, &nAtomType);
114 fscanf(fpSTATE, "%s %d", dummy, &nBondType);
115 fscanf(fpSTATE, "%s %lf", dummy, &region[1]);
116 fscanf(fpSTATE, "%s %lf", dummy, &region[2]);
117
118 if(timeNow == 0.0) region[2] *= 1.5; //Remove this when put on GitHub
119
120 density = nAtom/(region[1]*region[2]);
121 cells[1] = region[1] / rCut;
122 cells[2] = region[2] / rCut;
123 cellList = (int*)malloc((nAtom + cells[1] * cells[2] + 1) * sizeof(int));
124 regionH[1] = 0.5*region[1];
125 regionH[2] = 0.5*region[2];
126
127 //strain information
131 int n;
132
133 rx = (double*)malloc((nAtom + 1) * sizeof(double));
134 ry = (double*)malloc((nAtom + 1) * sizeof(double));
135 vx = (double*)malloc((nAtom + 1) * sizeof(double));
136 vy = (double*)malloc((nAtom + 1) * sizeof(double));
137 ax = (double*)malloc((nAtom + 1) * sizeof(double));
138 ay = (double*)malloc((nAtom + 1) * sizeof(double));
139 fx = (double*)malloc((nAtom + 1) * sizeof(double));
140 fy = (double*)malloc((nAtom + 1) * sizeof(double));
141 fax = (double*)malloc((nAtom + 1) * sizeof(double));
142 fay = (double*)malloc((nAtom + 1) * sizeof(double));
143 atomID = (int*)malloc((nAtom+1) * sizeof(int));
144 atomType = (int*)malloc((nAtom+1) * sizeof(int));
145 atomRadius = (double*)malloc((nAtom + 1) * sizeof(double));
146 atomMass = (double*)malloc((nAtom + 1) * sizeof(double));
147 speed = (double*)malloc((nAtom + 1) * sizeof(double));
148 discDragx = (double*)malloc((nAtom + 1) * sizeof(double));
149 discDragy = (double*)malloc((nAtom + 1) * sizeof(double));
150 molID = (int*)malloc((nAtom+1) * sizeof(int));
151
152 BondID = (int*)malloc((nBond+1)*sizeof(int));
153 BondType = (int*)malloc((nBond+1)*sizeof(int));
154 atom1 = (int*)malloc((nBond+1)*sizeof(int));
155 atom2 = (int*)malloc((nBond+1)*sizeof(int));
156 kb = (double*)malloc((nBond+1)*sizeof(double));
157 ro = (double*)malloc((nBond+1)*sizeof(double));
158 BondEnergy = (double*)malloc((nBond+1)*sizeof(double));
159 BondLength =(double*)malloc((nBond+1)*sizeof(double));
160 nodeDragx = (double*)malloc((nAtom + 1) * sizeof(double));
161 nodeDragy = (double*)malloc((nAtom + 1) * sizeof(double));
162 rxUnwrap = (double*)malloc((nAtom + 1) * sizeof(double));
163 ryUnwrap = (double*)malloc((nAtom + 1) * sizeof(double));
164 ImageX = (int*)malloc((nAtom+1) * sizeof(int));
165 ImageY = (int*)malloc((nAtom+1) * sizeof(int));
166
167
168 for(n = 1; n <= nAtom; n ++){
169 atomMass[n] = mass;
170 }
171
172 fscanf(fpSTATE, "%s\n", dummy);
173 for(n = 1; n <= nAtom; n ++)
174 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]);
175
176
177 fscanf(fpSTATE, "%s\n", dummy);
178 for(n=1; n<=nBond; n++)
179 fscanf(fpSTATE, "%d %d %d %d %lf %lf\n", &BondID[n], &BondType[n], &atom1[n], &atom2[n], &kb[n], &ro[n]);
180
181 fclose(fpSTATE);
182
183 //2D-List of bonded atoms. This is used to remove pair interaction
184 //calculation for the bonded atoms
185 isBonded = (int**)malloc((nAtom + 1) * sizeof(int*));
186 for (int i = 0; i <= nAtom; i++) {
187 isBonded[i] = (int*)malloc((nAtom + 1) * sizeof(int));
188 for (int j = 0; j <= nAtom; j++) {
189 isBonded[i][j] = 0;
190 }
191 }
192
193 for (n = 1; n <= nBond; n++) {
194 int i = atom1[n];
195 int j = atom2[n];
196 isBonded[i][j] = 1;
197 isBonded[j][i] = 1; // symmetric
198}
199
200 // List the interface atoms
201 nAtomInterface = 0;
202 nAtomBlock = 0;
203 nDiscInterface = 0;
204 bigDiameter = 2.8;
206
207 for(n = 1; n <= nAtom; n++){
208 if(fabs(ry[n]) < InterfaceWidth){
210 }
211 if(molID[n] == 2){
212 nAtomBlock++;
213 }
214 if(atomRadius[n] != 0.0){
216 } }
217
218
219 int BondPairInteract;
220 BondPairInteract = 0;
221 int m;
222 if(BondPairInteract == 1){
223 atomIDInterface = (int *)malloc((nAtomInterface+1)*sizeof(int));
224 m = 1;
225 for(n=1; n<=nAtom; n++){
226 if(fabs(ry[n]) < InterfaceWidth){
227 atomIDInterface[m] = atomID[n];
228 m++;
229 } } }
230 else if(BondPairInteract == 0){
232 atomIDInterface = (int *)malloc((nAtomInterface+1)*sizeof(int));
233 m = 1;
234 for(n=1; n<=nAtom; n++){
235 if(atomRadius[n] != 0.0){
236 atomIDInterface[m] = atomID[n];
237 m++;
238 } } }
239
241 PairID = (int*)malloc((nPairTotal+1) * sizeof(int));
242 Pairatom1 = (int*)malloc((nPairTotal+1) * sizeof(int));
243 Pairatom2 = (int*)malloc((nPairTotal+1) * sizeof(int));
244 PairXij = (double*)malloc((nPairTotal+1) * sizeof(double));
245 PairYij = (double*)malloc((nPairTotal+1) * sizeof(double));
246
247 fprintf(fpresult, "------------------------------------\n");
248 fprintf(fpresult, "-------PARAMETERS-----------\n");
249 fprintf(fpresult, "------------------------------------\n");
250 fprintf(fpresult, "nAtom\t\t\t%d\n", nAtom);
251 fprintf(fpresult, "nBond\t\t\t%d\n", nBond);
252 fprintf(fpresult, "nAtomType\t\t%d\n", nAtomType);
253 fprintf(fpresult, "nBondType\t\t%d\n", nBondType);
254 fprintf(fpresult, "nAtomBlock\t\t%d\n", nAtomBlock);
255 fprintf(fpresult, "nAtomInterface\t\t%d\n", nAtomInterface);
256 fprintf(fpresult, "nDiscInterface\t\t%d\n", nDiscInterface);
257 fprintf(fpresult, "mass\t\t\t%0.6g\n", mass);
258 fprintf(fpresult, "gamman\t\t\t%0.6g\n", gamman);
259 fprintf(fpresult, "strain\t\t\t%0.6g\n", strain);
260 fprintf(fpresult, "strainRate\t\t%0.6g\n", strainRate);
261 fprintf(fpresult, "FyBylx\t\t\t%0.6g\n", FyBylx);
262 fprintf(fpresult, "fxByfy\t\t\t%0.6g\n", fxByfy);
263 fprintf(fpresult, "DeltaY\t\t\t%0.6g\n", DeltaY);
264 fprintf(fpresult, "DeltaX\t\t\t%0.6g\n", DeltaX);
265 fprintf(fpresult, "HaltCondition\t\t%0.6g\n", HaltCondition);
266 fprintf(fpresult, "kappa\t\t\t%g\n", kappa);
267 fprintf(fpresult, "density\t\t\t%g\n", density);
268 fprintf(fpresult, "rCut\t\t\t%g\n", rCut);
269 fprintf(fpresult, "deltaT\t\t\t%g\n", deltaT);
270 fprintf(fpresult, "stepEquil\t\t%d\n", stepEquil);
271 fprintf(fpresult, "stepLimit\t\t%d\n", stepLimit);
272 fprintf(fpresult, "region[1]\t\t%0.16lf\n", region[1]);
273 fprintf(fpresult, "region[2]\t\t%0.16lf\n", region[2]);
274 fprintf(fpresult, "cells[1]\t\t%d\n", cells[1]);
275 fprintf(fpresult, "cells[2]\t\t%d\n", cells[2]);
276 fprintf(fpresult, "solver\t\t\t%s\n", solver);
277 fprintf(fpresult, "boundary\t\t%s %s\n", xBoundary, yBoundary);
278 fprintf(fpresult, "DampFlag\t\t%d\n", DampFlag);
279 fprintf(fpresult, "------------------------------------\n");
280 fprintf(fpresult, "#TimeNow TotalMomentum PotEngyPerAtom KinEngyPerAtom TotEngyPerAtom PairEnergyPerAtom BondEnergyPerAtom Press VirialSum\n");
281 fprintf(fpvrms, "#timeNow\tVrms \n");
282 fprintf(fpcom, "#timeNow\tComX\tComY\n");
283 fprintf(fpforce, "#timeNow\tforceSumxAtomType1\tforceSumyAtomType1\tforceSumxAtomType2\tforceSumyAtomType2\tforceSumxAtomType3\tforceSumyAtomType3\tforceSumxAtomType4\tforceSumyAtomType4\tforceSumxAtomType5\tforceSumyAtomType5\tforceSumxExtern\tforceSumyExtern\n");
284
285/* //Uncomment the following as per your acquirement
286 fprintf(fpstress, "strain %lf\n", strain);
287 fprintf(fpstress, "region[1] %lf\n", region[1]);
288 fprintf(fpstress, "region[2] %lf\n", region[2]);
289 fprintf(fpstress, "#timeNow virSumxx virSumyy virSumxy pressure\n");
290 fprintf(fpmomentum, "#timeNow Px Py\n");
291*/
292
293 if((strcmp(xBoundary, "p") != 0 && strcmp(xBoundary, "r") != 0) ||
294 (strcmp(yBoundary, "p") != 0 && strcmp(yBoundary, "r") != 0)) {
295 fprintf(fpresult, "Error: Invalid boundary value detected: '%s %s'. Only 'p' or 'r' are allowed.\n", xBoundary, yBoundary);
296 exit(EXIT_FAILURE); // Exit with failure status
297 }
298
299}
void Init()
Definition Init.c:31
void ReadBinaryRestart(const char *filename)
double kappa
Definition global.h:21
int limitCorrAv
Definition global.h:95
int nAtom
Definition global.h:24
int * BondID
double rCut
Definition global.h:21
double * discDragy
Definition global.h:42
int nDiscInterface
Definition global.h:76
double * atomMass
int nValCorr
Definition global.h:95
int cells[2+1]
Definition global.h:83
double FyBylx
Definition global.h:53
int nAtomBlock
Definition global.h:76
int nPairTotal
double mass
int * atom2
Definition global.h:36
int stepTraj
Definition global.h:24
double * fay
Definition global.h:85
int * BondType
Definition global.h:37
double DeltaY
int nBondType
Definition global.h:35
double Kn
int limitRdf
Definition global.h:104
double * speed
int nAtomType
double InterfaceWidth
double region[2+1]
double strain
int nValAcf
int stepRdf
Definition global.h:104
int stepCorr
Definition global.h:95
double gamman
int ** isBonded
int nBuffAcf
Definition global.h:100
int * molID
int * cellList
char yBoundary[10]
Definition global.h:68
double * ay
Definition global.h:17
int stepLimit
Definition global.h:24
double * vx
Definition global.h:17
int * ImageX
int * ImageY
Definition global.h:50
double density
Definition global.h:20
int stepAcf
Definition global.h:100
int nAtomInterface
double * ro
Definition global.h:38
double * BondLength
Definition global.h:39
int freezeAtomType
FILE * fpforce
FILE * fpresult
int * atomIDInterface
double deltaT
Definition global.h:20
double HaltCondition
int * atom1
double * rx
double * ax
Definition global.h:17
double * fy
Definition global.h:17
int nBond
double shearVelocity
Definition global.h:45
double * BondEnergy
double bigDiameter
Definition global.h:75
double * kb
double shearDisplacement
int stepEquil
Definition global.h:24
int stepAvg
Definition global.h:24
FILE * fpvrms
int DampFlag
int * PairID
int limitAcfAv
Definition global.h:100
double * PairXij
char mode[64]
FILE * fpcom
double * ryUnwrap
Definition global.h:51
double * nodeDragy
Definition global.h:42
int nBuffCorr
Definition global.h:95
double * atomRadius
char solver[128]
double timeNow
Definition global.h:20
int * Pairatom2
Definition global.h:63
char xBoundary[10]
int sizeHistRdf
Definition global.h:104
int stepCount
Definition global.h:24
char inputConfig[128]
Definition global.h:58
double rangeRdf
Definition global.h:103
double * PairYij
Definition global.h:64
double * discDragx
int * atomType
double strainRate
Definition global.h:44
double * vy
Definition global.h:17
double * rxUnwrap
int stepDump
Definition global.h:24
double fxByfy
Definition global.h:53
double regionH[2+1]
Definition global.h:20
double * nodeDragx
Definition global.h:42
int nFunCorr
Definition global.h:95
int * atomID
double DeltaX
Definition global.h:49
double * fx
Definition global.h:17
int * Pairatom1
Definition global.h:63
double * ry
Definition global.h:17
double * fax