-
Notifications
You must be signed in to change notification settings - Fork 0
/
mcmd.sct
executable file
·261 lines (260 loc) · 8.31 KB
/
mcmd.sct
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
#! /bin/bash
# dependencias:
CPUs=24
# RASPA, GULP
loc=`pwd`
objsMD=$loc/objsMD
objsMC=$loc/objsMC
#
mc_cycles=1000
md_cycles=52000
equi_cycles=48000
#
cation=$1
ncation=$2
temperature=$4
framework=$3
pressure=0.0
optimisate=yes
n_mc=10
n_md=5
RandomSeed='sun11jan13'
function print_input_RASPA {
# input on RASPA
let "movies_cycles = ncycles - 1"
echo "SimulationType $SimulationType
NumberOfCycles $ncycles
NumberOfInitializationCycles $inicialization
NumberOfEquilibrationCycles $equilibration
RestartFile $RestartFile
PrintEvery 100
Movies yes
WriteMoviesEvery $movies_cycles
RemoveAtomNumberCodeFromLabel yes
UseChargesFromCIFFile no
RemoveBondNeighboursFromLongRangeInteraction yes
RemoveBendNeighboursFromLongRangeInteraction yes
RemoveTorsionNeighboursFromLongRangeInteraction yes
InternalFrameworkLennardJonesInteractions no
Ensemble ${ensemble}
Forcefield local
CutOff 12.0
ChargeMethod Ewald
EwaldPrecision 1e-10
Framework 0
FrameworkName RHO_${framework}_$nAl
UnitCells 1 1 1
InputFileType cif
ExternalTemperature $temp
FlexibleFramework $FlexibleFramework
FrameworkDefinitions local
component 0 MoleculeName $cation
MoleculeDefinition local
ExtraFrameworkMolecule yes
TranslationProbability 1.0
RandomTranslationProbability 1.0
CreateNumberOfMolecules $ncations" > simulation.input
}
function inicialization_sim {
folder=ini
mkdir $folder
cd $folder
FrameworkName=RHO_${framework}_$nAl
cp $objsMC/*.def .
cp $objsMC/$FrameworkName.cif .
SimulationType=MC
RestartFile=no
ncations=$ncation
ncycles=0
inicialization=1
equilibration=1
FlexibleFramework=no
ensemble=NVT
temp=$temperature
print_input_RASPA
echo "inicialisation [RASPA]"
simulate
update
cd ..
}
function MC_move {
folder=MC_${i}_${j}
mkdir $folder
cd $folder
cp $objsMC/*.def .
cp $objsMD/* .
cp ../RHO_${framework}_$nAl.cif .
SimulationType=MC
RestartFile=no
FrameworkName=RHO_${framework}_$nAl
ncations=$ncation
ncycles=$mc_cycles
inicialization=1
equilibration=1
FlexibleFramework=yes
temp=$temperature
print_input_RASPA
echo "MC Cycle $i$j [RASPA]"
sed -i '/shO/d' RHO_${framework}_${nAl}.cif
simulate
sed '/MODEL 2/,$!d' Movies/System_0/Movie_*_allcomponents.pdb > c
sed 's/MODEL 2/MODEL 1/g' c > input.pdb
cp input.pdb ../input_${i}_${j}.pdb
rm c
gfortran -O2 -march=native pdb2gin.f90 -o pdb2gin
echo "Sorting PDB"
head -n2 input.pdb > pdb.top
tail -n1 input.pdb > pdb.bottom
touch pdb.middle
for k in 'Si' 'Al' 'Li' 'Na' 'K' 'Cs' 'Ca' 'Sr' 'Fe' 'O ' 'sh' ; do
grep "$k " input.pdb | grep 'ATOM' >> pdb.middle
sed -i "/$k /d" input.pdb
done
cat pdb.top > input.pdb
cat pdb.middle >> input.pdb
cat pdb.bottom >> input.pdb
echo "PDB 2 gin"
./pdb2gin > /dev/null
sed -i 's/sh core/ O shel/g' gin
sed -i 's/sh core/ O shel/g' gin
cat potentials.lib >> gin
echo "optimise shell [ GULP ]"
if [ $optimisate == "yes" ] ; then
cat gin.top.opti.shell > opti.gin
cat gin >> opti.gin
cat gin.bottom.opti >> opti.gin
mpirun --np $CPUs gulp < opti.gin > opti.gout
fi
sleep 1
cp gin.res ../gin.res.$i.$j
cd ..
}
function choose_ener {
gfortran -O2 -march=native pdb2gin.f90 -o pdb2gin
energia=6006919099999.9
confi=0
for j in `seq 1 $n_mc` ; do
ener=`grep 'Current total potential energy:' ../MC_${i}_${j}/Output/System_0/*.data | grep 'avg.' | tail -n1 | awk '{print $8}' | sed 's/)//g'`
echo $j $ener
if [ `echo "$ener < $energia" | bc -l` == 1 ] ; then
confi=$j
energia=$ener
fi
done
echo "Select: MC [ $confi ], [ $energia ]"
}
function choose_boltz {
gfortran -O2 -march=native pdb2gin.f90 -o pdb2gin
gfortran -O2 -march=native choose_structure_mc.f90 -o choose_structure_mc
config=0
rm choose_structure_mc.input toto
echo $n_mc > choose_structure_mc.input
echo $temp >> choose_structure_mc.input
for j in `seq 1 $n_mc` ; do
ener=`grep 'Current total potential energy:' ../MC_${i}_${j}/Output/System_0/*.data | grep 'avg.' | tail -n1 | awk '{print $8}' | sed 's/)//g'`
echo "$j $ener" >> toto
done
sort -nk2 toto >> choose_structure_mc.input
./choose_structure_mc < choose_structure_mc.input > choose_structure_mc.output
cat choose_structure_mc.output
confi=`tail -n1 choose_structure_mc.output | awk '{print $1}'`
energia=`tail -n1 choose_structure_mc.output | awk '{print $2}'`
echo "Select: MC [ $confi ], [ $energia ]"
}
function choose_boltz_gulp {
gfortran -O2 -march=native choose_structure_mc.f90 -o choose_structure_mc
config=0
rm choose_structure_mc.input toto
echo $n_mc > choose_structure_mc.input
echo $temp >> choose_structure_mc.input
for j in `seq 1 $n_mc` ; do
ener=`grep 'totalenergy' ../gin.res.$i.$j | awk '{print $2}'`
ener=`echo "$ener * 100000.0 / 8.6173324" | bc -l`
echo "$j $ener" >> toto
done
sort -nk2 toto >> choose_structure_mc.input
./choose_structure_mc < choose_structure_mc.input > choose_structure_mc.output
cat choose_structure_mc.output
confi=`tail -n1 choose_structure_mc.output | awk '{print $1}'`
energia=`tail -n1 choose_structure_mc.output | awk '{print $2}'`
echo "Select: MC [ $confi ], [ $energia ]"
}
#
function MD_move {
timestep=0.0002
temp=$temperature
ncycles=$md_cycles
production=`echo "$timestep * $ncycles" | bc -l`
equilibration=`echo "$timestep * $equi_cycles" | bc -l`
produequi=`echo "$production + $equilibration" | bc -l`
folder=MD_$i
mkdir $folder
cd $folder
cp $objsMD/* .
choose_boltz_gulp
cp ../gin.res.$i.$confi gin.res
finalizado=0
while [ `echo "$finalizado == 0" | bc -l` == 1 ] ; do
# molecular dynamics
sed 's/opti conv/md conp/g' gin.res > md.gin
sed -i '/output/d' md.gin
sed -i '/dump every/d' md.gin
cat gin.bottom.md >> md.gin
sed -i "s/TEMPERATURE/$temp/g" md.gin
sed -i "s/PRODUCTION/$production/g" md.gin
sed -i "s/EQUILIBRATION/$equilibration/g" md.gin
produequi=`echo "$production + $equilibration" | bc -l`
sed -i "s/PRODUEQUI/$produequi/g" md.gin
sed -i "s/TIMESTEP/$timestep/g" md.gin
mpirun --np $CPUs gulp < md.gin > md.gout.$i
finalizado=`grep 'Job Finished' md.gout.${i} | wc -l`
if [ `echo "$finalizado == 1" | bc -l` == 1 ] ; then
echo "Simulation finish"
touch GREAT
else
echo "Simulation wrong"
touch ERROR
confi=`echo $RANDOM % $n_mc + 1 | bc`
cp ../gin.res.$i.$j gin.res
cp ../input_${i}_${confi}.pdb input.pdb
fi
done
sed 's/md conp/single noenergy/g' md.res > res2cif.gin
gulp < res2cif.gin > res2cif.gout
# -----------------------------------
cp md.cif out.cif
cp out.cif out.cif.$i
for k in 'Li' 'Na' 'K' 'Cs' 'Ca' 'Sr' 'Fe' ; do sed -i "/$k /d" out.cif ; done
cp out.cif ../RHO_${framework}_${nAl}_$i.cif
cp out.cif ../RHO_${framework}_${nAl}.cif
cd ..
}
#
function update {
if [ -d RestartInitial ] ; then rm -R RestartInitial ; fi
cp -R Restart RestartInitial
cp Movies/System_0/Framework_0_final_1_1_1_P1.cif ../RHO_${framework}_${nAl}.cif
}
# MAIN: sodium, lithium, potassium, calcium, strontium
if [[ $cation = *sodium* ]] || [[ $cation = *lithium* ]] || [[ $cation == *potassium* ]] || [[ $cation == *caesium* ]] ; then
nAl=$ncation
elif [[ $cation = *calcium* ]] || [[ $cation = *strontium* ]] || [[ $cation == *magnesium* ]] ; then
nAl=`echo "$ncation * 2.0" | bc -l | sed -e 's/\([0-9]\+\)\.[0-9]\+/\1/g'`
else
nAl=$ncation
fi
folder=${cation}_${temperature}_${RandomSeed}_${framework}_$nAl
if [ ! -d $folder ] ; then
mkdir $folder
cd $folder
inicialization_sim
for i in `seq 1 $n_md` ; do
for j in `seq 1 $n_mc` ; do
MC_move
done
temp=$temperature
MD_move
done
cd ..
fi
exit 0