Lemina
A molecular dynamics package for network, granular material and point particles with a range of interaction potential.
 
Loading...
Searching...
No Matches
ReadRestartBinary.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#include <stdio.h>
22#include <stdlib.h>
23#include <string.h>
24#include"global.h"
25
26// Must match RestartHeader definition in WriteBinaryRestart
27typedef struct {
28 char magic[8]; // "LAMINA\0"
29 double version;
30 double timeNow;
31 int nAtom;
32 int nBond;
35 double regionX; // = region[1]
36 double regionY; // = region[2]
44 double uSumPair;
45 double virSumPair;
50 double virSumBond;
58
59void ReadBinaryRestart(const char *filename) {
60
61 FILE *fp = fopen(filename, "rb");
62 if (!fp) {
63 fprintf(stderr, "Error opening binary restart file %s for reading\n", filename);
64 exit(EXIT_FAILURE);
65 }
66 RestartHeader hdr; //Declare here
67 fread(&hdr, sizeof(RestartHeader), 1, fp); //Use it
68
69 if(strncmp(hdr.magic, "LAMINA", 6) != 0) {
70 fprintf(stderr, "Invalid file format: magic = %.8s [from %s()]\n", hdr.magic, __func__);
71 fclose(fp);
72 exit(EXIT_FAILURE); //Must return void, not return 1
73 }
74
75 //Now assigned the values that were read from binary file to global parameters
76 timeNow = hdr.timeNow;
77 nAtom = hdr.nAtom;
78 nBond = hdr.nBond;
79 nAtomType = hdr.nAtomType;
80 nBondType = hdr.nBondType;
81 region[1] = hdr.regionX;
82 region[2] = hdr.regionY;
90 uSumPair = hdr.uSumPair;
100 stepCount = hdr.stepCount;
103
104 density = nAtom / (region[1] * region[2]);
105 cells[1] = region[1] / rCut;
106 cells[2] = region[2] / rCut;
107 regionH[1] = 0.5 * region[1];
108 regionH[2] = 0.5 * region[2];
109
110 cellList = (int *)malloc((nAtom + cells[1] * cells[2] + 1) * sizeof(int));
111
112 printf("Running from restart file:\n");
113 printf("Running: %s, version: %0.3lf\n", hdr.magic, hdr.version);
114 printf("timeNow: %lf\n", timeNow);
115 printf("stepCount: %d\n", stepCount);
116
117 //Allocating the memory to arrays
118 atomID = (int *)malloc((nAtom + 1) * sizeof(int));
119 molID = (int *)malloc((nAtom + 1) * sizeof(int));
120 atomType = (int *)malloc((nAtom + 1) * sizeof(int));
121 atomRadius = (double *)malloc((nAtom + 1) * sizeof(double));
122 rx = (double *)malloc((nAtom + 1) * sizeof(double));
123 ry = (double *)malloc((nAtom + 1) * sizeof(double));
124 vx = (double *)malloc((nAtom + 1) * sizeof(double));
125 vy = (double *)malloc((nAtom + 1) * sizeof(double));
126 ax = (double *)malloc((nAtom + 1) * sizeof(double));
127 ay = (double *)malloc((nAtom + 1) * sizeof(double));
128 fx = (double *)malloc((nAtom + 1) * sizeof(double));
129 fy = (double *)malloc((nAtom + 1) * sizeof(double));
130 atomMass = (double *)malloc((nAtom + 1) * sizeof(double));
131 discDragx = (double *)malloc((nAtom + 1) * sizeof(double));
132 discDragy = (double *)malloc((nAtom + 1) * sizeof(double));
133 atomIDInterface = (int *)malloc((nAtom + 1) * sizeof(int));
134
135 BondID = (int *)malloc((nBond + 1) * sizeof(int));
136 BondType = (int *)malloc((nBond + 1) * sizeof(int));
137 atom1 = (int *)malloc((nBond + 1) * sizeof(int));
138 atom2 = (int *)malloc((nBond + 1) * sizeof(int));
139 kb = (double *)malloc((nBond + 1) * sizeof(double));
140 ro = (double *)malloc((nBond + 1) * sizeof(double));
141 BondEnergy = (double *)malloc((nBond + 1) * sizeof(double));
142 BondLength = (double *)malloc((nBond + 1) * sizeof(double));
143 nodeDragx = (double *)malloc((nAtom + 1) * sizeof(double));
144 nodeDragy = (double *)malloc((nAtom + 1) * sizeof(double));
145 rxUnwrap = (double *)malloc((nAtom + 1) * sizeof(double));
146 ryUnwrap = (double *)malloc((nAtom + 1) * sizeof(double));
147 ImageX = (int *)malloc((nAtom + 1) * sizeof(int));
148 ImageY = (int *)malloc((nAtom + 1) * sizeof(int));
149
150 PairID = (int *)malloc((nPairTotal + 1) * sizeof(int));
151 Pairatom1 = (int *)malloc((nPairTotal + 1) * sizeof(int));
152 Pairatom2 = (int *)malloc((nPairTotal + 1) * sizeof(int));
153 PairXij = (double *)malloc((nPairTotal + 1) * sizeof(double));
154 PairYij = (double *)malloc((nPairTotal + 1) * sizeof(double));
155
156 //Here we are reading the data to binary file
157 fread(&atomID[1], sizeof(int), nAtom, fp);
158 fread(&molID[1], sizeof(int), nAtom, fp);
159 fread(&atomType[1], sizeof(int), nAtom, fp);
160 fread(&atomRadius[1], sizeof(double), nAtom, fp);
161 fread(&rx[1], sizeof(double), nAtom, fp);
162 fread(&ry[1], sizeof(double), nAtom, fp);
163 fread(&vx[1], sizeof(double), nAtom, fp);
164 fread(&vy[1], sizeof(double), nAtom, fp);
165 fread(&ax[1], sizeof(double), nAtom, fp);
166 fread(&ay[1], sizeof(double), nAtom, fp);
167 fread(&fx[1], sizeof(double), nAtom, fp);
168 fread(&fy[1], sizeof(double), nAtom, fp);
169 fread(&atomMass[1], sizeof(double), nAtom, fp);
170 fread(&discDragx[1], sizeof(double), nAtom, fp);
171 fread(&discDragy[1], sizeof(double), nAtom, fp);
172 fread(&atomIDInterface[1], sizeof(int), nAtomInterface, fp);
173
174 fread(&BondID[1], sizeof(int), nBond, fp);
175 fread(&BondType[1], sizeof(int), nBond, fp);
176 fread(&atom1[1], sizeof(int), nBond, fp);
177 fread(&atom2[1], sizeof(int), nBond, fp);
178 fread(&kb[1], sizeof(double), nBond, fp);
179 fread(&ro[1], sizeof(double), nBond, fp);
180 fread(&BondEnergy[1], sizeof(double), nBond, fp);
181 fread(&BondLength[1], sizeof(double), nBond, fp);
182 fread(&nodeDragx[1], sizeof(double), nAtom, fp);
183 fread(&nodeDragy[1], sizeof(double), nAtom, fp);
184 fread(&rxUnwrap[1], sizeof(double), nAtom, fp);
185 fread(&ryUnwrap[1], sizeof(double), nAtom, fp);
186 fread(&ImageX[1], sizeof(int), nAtom, fp);
187 fread(&ImageY[1], sizeof(int), nAtom, fp);
188
189 fread(&PairID[1], sizeof(int), nPairActive, fp);
190 fread(&Pairatom1[1], sizeof(int), nPairActive, fp);
191 fread(&Pairatom2[1], sizeof(int), nPairActive, fp);
192 fread(&PairXij[1], sizeof(double), nPairActive, fp);
193 fread(&PairYij[1], sizeof(double), nPairActive, fp);
194
195 //2D-List of bonded atoms. This is used to remove pair interaction
196 //calculation for the bonded atoms
197 isBonded = (int **)malloc((nAtom + 1) * sizeof(int*));
198 for (int i = 0; i <= nAtom; i++) {
199 isBonded[i] = (int *)malloc((nAtom + 1) * sizeof(int));
200 for (int j = 0; j <= nAtom; j++) {
201 isBonded[i][j] = 0;
202 } }
203 for (int n = 1; n <= nBond; n++) {
204 int i = atom1[n];
205 int j = atom2[n];
206 isBonded[i][j] = 1;
207 isBonded[j][i] = 1; // symmetric
208}
209
210 fprintf(fpresult, "------------------------------------\n");
211 fprintf(fpresult, "-------PARAMETERS-----------\n");
212 fprintf(fpresult, "------------------------------------\n");
213 fprintf(fpresult, "nAtom\t\t\t%d\n", nAtom);
214 fprintf(fpresult, "nBond\t\t\t%d\n", nBond);
215 fprintf(fpresult, "nAtomType\t\t%d\n", nAtomType);
216 fprintf(fpresult, "nBondType\t\t%d\n", nBondType);
217 fprintf(fpresult, "nAtomBlock\t\t%d\n", nAtomBlock);
218 fprintf(fpresult, "nAtomInterface\t\t%d\n", nAtomInterface);
219 fprintf(fpresult, "nDiscInterface\t\t%d\n", nDiscInterface);
220 fprintf(fpresult, "mass\t\t\t%0.6g\n", mass);
221 fprintf(fpresult, "gamman\t\t\t%0.6g\n", gamman);
222 fprintf(fpresult, "strain\t\t\t%0.6g\n", strain);
223 fprintf(fpresult, "strainRate\t\t%0.6g\n", strainRate);
224 fprintf(fpresult, "FyBylx\t\t\t%0.6g\n", FyBylx);
225 fprintf(fpresult, "fxByfy\t\t\t%0.6g\n", fxByfy);
226 fprintf(fpresult, "DeltaY\t\t\t%0.6g\n", DeltaY);
227 fprintf(fpresult, "DeltaX\t\t\t%0.6g\n", DeltaX);
228 fprintf(fpresult, "HaltCondition\t\t%0.6g\n", HaltCondition);
229 fprintf(fpresult, "kappa\t\t\t%g\n", kappa);
230 fprintf(fpresult, "density\t\t\t%g\n", density);
231 fprintf(fpresult, "rCut\t\t\t%g\n", rCut);
232 fprintf(fpresult, "deltaT\t\t\t%g\n", deltaT);
233 fprintf(fpresult, "stepEquil\t\t%d\n", stepEquil);
234 fprintf(fpresult, "stepLimit\t\t%d\n", stepLimit);
235 fprintf(fpresult, "region[1]\t\t%0.16lf\n", region[1]);
236 fprintf(fpresult, "region[2]\t\t%0.16lf\n", region[2]);
237 fprintf(fpresult, "cells[1]\t\t%d\n", cells[1]);
238 fprintf(fpresult, "cells[2]\t\t%d\n", cells[2]);
239 fprintf(fpresult, "solver\t\t\t%s\n", solver);
240 fprintf(fpresult, "boundary\t\t%s %s\n", xBoundary, yBoundary);
241 fprintf(fpresult, "DampFlag\t\t%d\n", DampFlag);
242 fprintf(fpresult, "------------------------------------\n");
243 fprintf(fpresult, "#TimeNow TotalMomentum PotEngyPerAtom KinEngyPerAtom TotEngyPerAtom PairEnergyPerAtom BondEnergyPerAtom Press VirialSum\n");
244 fprintf(fpvrms, "#timeNow\tVrms \n");
245 fprintf(fpcom, "#timeNow\tComX\tComY\n");
246 fprintf(fpforce, "#timeNow\tforceSumxAtomType1\tforceSumyAtomType1\tforceSumxAtomType2\tforceSumyAtomType2\tforceSumxAtomType3\tforceSumyAtomType3\tforceSumxAtomType4\tforceSumyAtomType4\tforceSumxAtomType5\tforceSumyAtomType5\tforceSumxExtern\tforceSumyExtern\n");
247 fclose(fp);
248 }
249
void ReadBinaryRestart(const char *filename)
double kappa
Definition global.h:21
double virSumBond
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 cells[2+1]
Definition global.h:83
double FyBylx
Definition global.h:53
double virSumPairyy
Definition global.h:88
int nAtomBlock
Definition global.h:76
double virSumPair
Definition global.h:88
int nPairTotal
double mass
int * atom2
Definition global.h:36
int * BondType
Definition global.h:37
double DeltaY
int nBondType
Definition global.h:35
int nAtomType
double InterfaceWidth
double region[2+1]
double strain
double virSumPairxx
Definition global.h:88
int nPairActive
Definition global.h:62
double gamman
int ** isBonded
int * molID
double virSumBondxy
Definition global.h:89
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 nAtomInterface
double * ro
Definition global.h:38
double * BondLength
Definition global.h:39
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 * BondEnergy
double bigDiameter
Definition global.h:75
double * kb
int stepEquil
Definition global.h:24
FILE * fpvrms
int DampFlag
double uSumPair
int * PairID
double * PairXij
FILE * fpcom
double * ryUnwrap
Definition global.h:51
double * nodeDragy
Definition global.h:42
double * atomRadius
double virSumBondxx
Definition global.h:89
char solver[128]
double timeNow
Definition global.h:20
int * Pairatom2
Definition global.h:63
char xBoundary[10]
int stepCount
Definition global.h:24
double forceSumyExtern
Definition global.h:53
double TotalBondEnergy
double * PairYij
Definition global.h:64
double * discDragx
int * atomType
double strainRate
Definition global.h:44
double * vy
Definition global.h:17
double * rxUnwrap
double fxByfy
Definition global.h:53
double regionH[2+1]
Definition global.h:20
double * nodeDragx
Definition global.h:42
double virSumBondyy
Definition global.h:89
int * atomID
double forceSumxExtern
Definition global.h:53
double DeltaX
Definition global.h:49
double virSumPairxy
Definition global.h:88
double * fx
Definition global.h:17
int * Pairatom1
Definition global.h:63
double * ry
Definition global.h:17