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
29
void
ReadBinaryRestart
(
const
char
*filename);
30
31
void
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
128
strainRate
=
strain
/
deltaT
;
129
shearDisplacement
=
strain
*
region
[2];
130
shearVelocity
=
strainRate
*
region
[2];
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;
205
InterfaceWidth
= 5.0 *
bigDiameter
;
206
207
for
(n = 1; n <=
nAtom
; n++){
208
if
(fabs(
ry
[n]) <
InterfaceWidth
){
209
nAtomInterface
++;
210
}
211
if
(
molID
[n] == 2){
212
nAtomBlock
++;
213
}
214
if
(
atomRadius
[n] != 0.0){
215
nDiscInterface
++;
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){
231
nAtomInterface
=
nDiscInterface
;
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
240
nPairTotal
= 0.5 *
nAtomInterface
* (
nAtomInterface
-1);
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
}
Init
void Init()
Definition
Init.c:31
ReadBinaryRestart
void ReadBinaryRestart(const char *filename)
Definition
ReadRestartBinary.c:59
global.h
kappa
double kappa
Definition
global.h:21
limitCorrAv
int limitCorrAv
Definition
global.h:95
nAtom
int nAtom
Definition
global.h:24
BondID
int * BondID
rCut
double rCut
Definition
global.h:21
discDragy
double * discDragy
Definition
global.h:42
nDiscInterface
int nDiscInterface
Definition
global.h:76
atomMass
double * atomMass
nValCorr
int nValCorr
Definition
global.h:95
cells
int cells[2+1]
Definition
global.h:83
FyBylx
double FyBylx
Definition
global.h:53
nAtomBlock
int nAtomBlock
Definition
global.h:76
nPairTotal
int nPairTotal
mass
double mass
atom2
int * atom2
Definition
global.h:36
stepTraj
int stepTraj
Definition
global.h:24
fay
double * fay
Definition
global.h:85
BondType
int * BondType
Definition
global.h:37
DeltaY
double DeltaY
nBondType
int nBondType
Definition
global.h:35
Kn
double Kn
limitRdf
int limitRdf
Definition
global.h:104
speed
double * speed
nAtomType
int nAtomType
InterfaceWidth
double InterfaceWidth
region
double region[2+1]
strain
double strain
nValAcf
int nValAcf
stepRdf
int stepRdf
Definition
global.h:104
stepCorr
int stepCorr
Definition
global.h:95
gamman
double gamman
isBonded
int ** isBonded
nBuffAcf
int nBuffAcf
Definition
global.h:100
molID
int * molID
cellList
int * cellList
yBoundary
char yBoundary[10]
Definition
global.h:68
ay
double * ay
Definition
global.h:17
stepLimit
int stepLimit
Definition
global.h:24
vx
double * vx
Definition
global.h:17
ImageX
int * ImageX
ImageY
int * ImageY
Definition
global.h:50
density
double density
Definition
global.h:20
stepAcf
int stepAcf
Definition
global.h:100
nAtomInterface
int nAtomInterface
ro
double * ro
Definition
global.h:38
BondLength
double * BondLength
Definition
global.h:39
freezeAtomType
int freezeAtomType
fpforce
FILE * fpforce
fpresult
FILE * fpresult
atomIDInterface
int * atomIDInterface
deltaT
double deltaT
Definition
global.h:20
HaltCondition
double HaltCondition
atom1
int * atom1
rx
double * rx
ax
double * ax
Definition
global.h:17
fy
double * fy
Definition
global.h:17
nBond
int nBond
shearVelocity
double shearVelocity
Definition
global.h:45
BondEnergy
double * BondEnergy
bigDiameter
double bigDiameter
Definition
global.h:75
kb
double * kb
shearDisplacement
double shearDisplacement
stepEquil
int stepEquil
Definition
global.h:24
stepAvg
int stepAvg
Definition
global.h:24
fpvrms
FILE * fpvrms
DampFlag
int DampFlag
PairID
int * PairID
limitAcfAv
int limitAcfAv
Definition
global.h:100
PairXij
double * PairXij
mode
char mode[64]
fpcom
FILE * fpcom
ryUnwrap
double * ryUnwrap
Definition
global.h:51
nodeDragy
double * nodeDragy
Definition
global.h:42
nBuffCorr
int nBuffCorr
Definition
global.h:95
atomRadius
double * atomRadius
solver
char solver[128]
timeNow
double timeNow
Definition
global.h:20
Pairatom2
int * Pairatom2
Definition
global.h:63
xBoundary
char xBoundary[10]
sizeHistRdf
int sizeHistRdf
Definition
global.h:104
stepCount
int stepCount
Definition
global.h:24
inputConfig
char inputConfig[128]
Definition
global.h:58
rangeRdf
double rangeRdf
Definition
global.h:103
PairYij
double * PairYij
Definition
global.h:64
discDragx
double * discDragx
atomType
int * atomType
strainRate
double strainRate
Definition
global.h:44
vy
double * vy
Definition
global.h:17
rxUnwrap
double * rxUnwrap
stepDump
int stepDump
Definition
global.h:24
fxByfy
double fxByfy
Definition
global.h:53
regionH
double regionH[2+1]
Definition
global.h:20
nodeDragx
double * nodeDragx
Definition
global.h:42
nFunCorr
int nFunCorr
Definition
global.h:95
atomID
int * atomID
DeltaX
double DeltaX
Definition
global.h:49
fx
double * fx
Definition
global.h:17
Pairatom1
int * Pairatom1
Definition
global.h:63
ry
double * ry
Definition
global.h:17
fax
double * fax
source
Init.c
Generated by
1.13.2