-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathshake.m
138 lines (112 loc) · 3.82 KB
/
shake.m
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
function positions = shake(lattice, parameters)
%SHAKE Randomizes positons of the particles.
% shake(lattice) uses Monte Carlo method to randomize positions of the
% sphere which initially were arranged in some crystalline structure.
% Initial length of neighbor list. It will be automatically extended if
% necesseery, so I don't see a point to make it an adjustable parameter for
% now.
nneighs = 12;
dmax = parameters.dmax;
ratiomin = parameters.acceptmin;
ratiomax = parameters.acceptmax;
npart = lattice.npart;
part.r = [0; 0; 0];
part.neighbors = -ones(1, nneighs);
particles = repmat(part, 1, npart);
for i = 1:npart
particles(i).r = lattice.positions(:, i);
end
box = lattice.box;
% Initialize neighbor list for each particle.
particles = findneigh(particles, box, parameters.rcut);
acceptacc = 0.0;
totalacc = 0.0;
ratio = 0.0;
ncycles = 1;
while (ratio < ratiomin || ratio > ratiomax || ncycles == parameters.maxcycles)
[particles, accepted] = makecycle(particles, box, dmax);
acceptacc = acceptacc + accepted;
totalacc = totalacc + npart;
% Adujst maximal displcement if the acceptance ratio is out of bounds.
if mod(ncycles, parameters.nadjust) == 0
ratio = acceptacc / totalacc;
if ratio < parameters.acceptmin
dmax = 0.99 * dmax;
end
if ratio > parameters.acceptmax
dmax = 1.01 * dmax;
end
% Reset counters.
acceptacc = 0.0;
totalacc = 0.0;
end
% Refresh neihbor lists from time to time.
if mod(ncycles, parameters.nupdate) == 0
for i = 1:npart
particles(i).neighbors = -ones(1, nneighs);
end
particles = findneigh(particles, box, parameters.rcut);
end
ncycles = ncycles + 1;
end
[ratio, dmax]
res = 100;
rdfvalues = zeros(1, res);
ndump = 0;
acceptacc = 0.0;
totalacc = 0.0;
for i = 1:parameters.maxcycles
[particles, accepted] = makecycle(particles, box, dmax);
acceptacc = acceptacc + accepted;
totalacc = totalacc + npart;
% Refresh neighbor lists from time to time.
if mod(i, parameters.nupdate) == 0
for j = 1:npart
particles(j).neighbors = -ones(1, nneighs);
end
particles = findneigh(particles, box, parameters.rcut);
end
% Collect data for RDF.
if mod(i, parameters.ndump) == 0
[values, centers] = rdf(particles, box, res);
rdfvalues = rdfvalues + values;
ndump = ndump + 1;
end
end
rdfcenters = centers;
rdfvalues = rdfvalues / (npart * ndump);
tab = [rdfcenters', rdfvalues'];
save('rdf.dat', 'tab', '-ascii');
positions = zeros(3, npart);
for i = 1:npart
positions(:, i) = particles(i).r;
end
end
function [particles, accepted] = makecycle(particles, box, maxdisp)
%MAKECYCLE Attempts to randomly move each particle from its position.
% [particles, accepted] = makecycle(particles, maxdisp) attempts to
% randomly generate a new position for each particle within a cube with
% edge length maxdisp centered on its current position.
%
% It returns particles with updated positions and number of successful
% attempts.
accepted = 0.0;
for i = 1:length(particles)
p = particles(i);
% Generate new tentative position.
p.r = p.r + maxdisp * (1.0 - 2.0 * rand(3, 1));
% Apply periodic boundary conditions along x, y (but NOT z) directions.
%p.r = p.r - [box(1); box(2); 0.0] .* fix(2.0 * p.r ./ box);
p.r = p.r - box .* fix(2.0 * p.r ./ box);
% Accept new position if it doesn't lead to overlaps with its
% neighbors and walls.
neighbors = particles(p.neighbors(p.neighbors > 0));
if isempty(neighbors)
warning('Particle %d has no neighbors.', i);
end
if ~isoverlap(p, neighbors, box)
particles(i) = p;
accepted = accepted + 1;
end
end
end