-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmain.cpp
211 lines (188 loc) · 7.4 KB
/
main.cpp
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
/****************************************************
** dislocMMC
**
** A Metropolis Monte Carlo algorithm to model the
** evolution of structure in defective sp2 bonded
** carbon systems
**
** The code calls the LAMMPS atomistic simulation
** program for structural optimisation with a
** reactive force-field
**
** T.Trevethan 2016
*****************************************************/
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <vector>
#include <random>
#include "bond.h"
#include "config.h"
#include "atom.h"
#include "polygon.h"
#include "setup.h"
int main()
{
std::cout << "dislocMMC v.0.60" << std::endl;
//initial configuration file
std::string inconfig = "data.in";
//options input file
std::string inoptions = "setup.in";
//name output files
std::string polyfile = "poly.xyz"; // 5-fold and 7-fold polygon coordinates: xyz movie
std::string outfile = "traj.xyz"; // full system coordinates: xyz movie
std::string bondfile = "bonds.xyz"; // active bond site coordinates: xyz movie
std::string enfile = "energy.out"; // optimised system energy per step
// polygon and active bond storage vectors
std::vector<Polygon> polys;
std::vector<Bond> bonds;
std::vector<Bond> bondc;
// get the simulation set-up and options and LAMMPS command from the input file
Setup setup(inoptions);
//initialise Mersenne Twister with seed
std::mt19937 generator (setup.seed());
// construct initial configuration from input file
Config config(inconfig);
// perform initial structure analysis
std::cout << "Initial structure analysis ..." << std::endl;
long numbnd = config.analyse(setup.coff(),bonds,bondc);
std::cout << "Found " << numbnd << " active bonds" << std::endl;
//write initial configuration and features (rotatable bonds and polygons) to output
config.write(outfile);
config.writePoly(polyfile);
config.writeBonds(bondfile);
//perform initial relaxation
std::cout << "First optimisation ..." << std::endl;
double en = config.relax(setup.lmp());
std::cout << "Complete. Final energy: " << en << " eV" << std::endl;
//do downhill energy minimisation
int ibnd = -1;
if(setup.type() == 1) {
for(int istep = 0; istep < setup.maxs(); istep++) {
//record each bond rotated configuration
std::vector<Config> configR;
configR.clear();
//loop over all rotatable atom pairs - optimise and find lowest energy
double minen = 0.0;
ibnd = -1;
for(int nb = 0; nb < numbnd; nb++) {
//copy configuration
Config tconfig = config;
//rotate bond
tconfig.rotate(nb);
std::cout << "Optimise configuration " << nb << " ..." << std::endl;
en = tconfig.relax(setup.lmp());
std::cout << "Complete. Final energy: " << en << " eV" << std::endl;
configR.push_back(tconfig);
//get the lowest energy structure
if(en < minen) {
minen = en;
ibnd = nb;
}
}
std::cout << "Selected rotation: " << ibnd << std::endl;
//copy lowest energy rotation into config
config = configR[ibnd];
numbnd = config.update(setup.coff(),ibnd);
std::cout << "Update. Bonds: " << numbnd << std::endl;
//write new configuration and features to output
config.write(outfile);
// config.writePoly(polyfile);
config.writeBonds(bondfile);
config.writeEn(enfile);
//update bond table
}
} else if(setup.type() == 2)
//do Metropolis-Hastings Monte Carlo
{
for(int istep = 0; istep < setup.maxs(); istep++) {
//pick a rotatable bond at random
double rnum = generator()*1.0/(generator.max()*1.0);
int nb = 0;
for(int ib = 0; ib < numbnd; ib++) {
double rscale = rnum*numbnd;
if(rscale > ib*1.0 && rscale <= (ib+1)*1.0) {
std::cout << "r " << rnum << " rs " << rscale << " ib " << ib << std::endl;
nb = ib;
break;
}
}
//copy configuration
Config tconfig = config;
//rotate bond
tconfig.rotate(nb);
std::cout << "Optimise configuration " << nb << " ..." << std::endl;
double en2 = tconfig.relax(setup.lmp());
std::cout << "Complete. Final energy: " << en << " eV" << std::endl;
//if energy is lower, accept move
if(en2 <= en) {
std::cout << "Energy change: " << (en2 - en) << " eV. Move accepted." << std::endl;
config = tconfig;
numbnd = config.update(setup.coff(),nb);
std::cout << "Update. Bonds: " << numbnd << std::endl;
//write new configuration and features to output
config.write(outfile);
// config.writePoly(polyfile);
config.writeBonds(bondfile);
config.writeEn(enfile);
} else {
//if energy is higher, calculate the probability from the Boltzman distribution
double sienergy = (en2 - en)*1.602e-19;
double boltz = exp(-sienergy/(1.38e-23*setup.tmp()));
std::cout << "Energy change: " << (en2 - en) << " eV. Boltz Prob: " << boltz << std::endl;
//get new random number and determine move
double rnum = generator()*1.0/(generator.max()*1.0);
if(boltz > rnum) {
std::cout << "Move Accepted." << std::endl;
config = tconfig;
numbnd = config.update(setup.coff(),nb);
std::cout << "Update. Bonds: " << numbnd << std::endl;
config.write(outfile);
// config.writePoly(polyfile);
config.writeBonds(bondfile);
config.writeEn(enfile);
} else {
std::cout << "Move Rejected." << std::endl;
}
}
}
} else if(setup.type() == 3)
//do double climb
{
config.dclimb(0);
config.relax(setup.lmp());
config.write(outfile);
} else if(setup.type() == 4)
//do single climb
{
config.sclimb(2);
config.relax(setup.lmp());
config.write(outfile);
config.analyse(setup.coff(),bonds,bondc);
std::vector<Polygon> octs;
long n8 = config.getOcts(octs);
std::cout << " 8s: " << n8 << std::endl;
config.shuffle(0,-1);
config.write(outfile);
config.relax(setup.lmp());
config.write(outfile);
config.shclimb(0);
config.relax(setup.lmp());
config.write(outfile);
} else if(setup.type() == 5)
//do single climb on shuffle
{
config.shclimb(0);
config.relax(setup.lmp());
config.write(outfile);
config.analyse(setup.coff(),bonds,bondc);
std::vector<Polygon> octs;
long n8 = config.getOcts(octs);
std::cout << " 8s: " << n8 << std::endl;
} else {
std::cout << "Invalid type option" << std::endl;
exit(1);
}
std::cout << "Complete. Exiting ..." << std::endl;
return 0;
}