-
Notifications
You must be signed in to change notification settings - Fork 0
/
mc_clsmc.cc
77 lines (66 loc) · 1.76 KB
/
mc_clsmc.cc
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
#include "mc_input.h"
#include "mc_confg.h"
#include "mc_setup.h"
#include "mc_utils.h"
#include "mc_poten.h"
#include "mc_clsmc.h"
#include "mc_randg.h"
#include <math.h>
// #define coodstest 1
double MCTotal;
double MCAccep;
void MCMove(void)
{
int numb=NumbMoveAtoms;
double disp[NDIM];
for(int atom=0;atom<numb;atom++)
{
int gatom=IMoving[atom];
#ifdef coodstest
cout<<"Moving : "<<gatom<<endl;
#endif
for (int id=0;id<NDIM;id++) // MOVE
{
disp[id] = MCAtom.mcstep*(rnd1()-0.5);
#ifdef coodstest
cout<<"Offset : "<<disp[id]<<endl;
#endif
newcoords[id][gatom] = MCCoords[id][gatom];
newcoords[id][gatom] += disp[id];
#ifdef coodstest
cout<<MCCoords[id][gatom]<<BLANK<<newcoords[id][gatom]<<endl;
#endif
}
double deltav = 0.0; // ACCEPT/REJECT
deltav += (MCPot(gatom,newcoords)-MCPot(gatom,MCCoords));
bool Accepted = false;
if (deltav<0.0) Accepted = true;
else if
(exp(-deltav/Temperature)>rnd2()) Accepted = true;
MCTotal+= 1.0;
if (Accepted)
{
MCAccep += 1.0;
for (int id=0;id<NDIM;id++) // save accepted configuration
MCCoords[id][gatom] = newcoords[id][gatom];
}
} // END sum over atoms (fixed atom type)
}
double MCPot(int atom0, double **pos)
{
double dr[NDIM];
double spot = 0.0;
for (int atom1=0;atom1<NumbAtoms;atom1++)
if (atom1 != atom0) // skip "self-interaction"
{
double dr2 = 0.0;
for (int id=0;id<NDIM;id++)
{
dr[id] = (pos[id][atom0] - MCCoords[id][atom1]);
dr2 += (dr[id]*dr[id]);
}
double r = sqrt(dr2);
spot += SPot1D(r); // 1D interaction
}
return (spot);
}