import colored langevin thermostat package from michele ceriotti

this is the unmodified patch without testing its functionality.
This commit is contained in:
Axel Kohlmeyer
2014-06-23 18:15:03 -04:00
parent 2e20b74d7c
commit 86d4432ff8
9 changed files with 2290 additions and 1 deletions

View File

@ -0,0 +1,412 @@
LAMMPS Description
192 atoms
128 bonds
64 angles
2 atom types
1 bond types
1 angle types
0 12.4244 xlo xhi
0 12.4244 ylo yhi
0 12.4244 zlo zhi
Masses
1 15.9994
2 1.0080
Atoms
1 1 1 -1.1128 2.7475 12.8896 2.6639
2 1 2 0.5564 3.1049 12.7379 1.8224
3 1 2 0.5564 1.7943 13.1205 2.5592
4 2 1 -1.1128 -0.7565 -1.1504 -1.5704
5 2 2 0.5564 -1.0109 -0.3411 -2.0812
6 2 2 0.5564 -0.7469 -1.7213 -2.2771
7 3 1 -1.1128 5.7663 -0.4112 0.4047
8 3 2 0.5564 6.4944 -0.7495 -0.3433
9 3 2 0.5564 4.9901 -0.2595 -0.2404
10 4 1 -1.1128 12.5178 13.4586 2.5833
11 4 2 0.5564 11.9595 12.7227 2.5625
12 4 2 0.5564 12.3870 13.9685 1.6020
13 5 1 -1.1128 11.5658 1.2126 -3.4671
14 5 2 0.5564 10.7292 1.5210 -4.0420
15 5 2 0.5564 12.2784 0.9456 -4.0624
16 6 1 -1.1128 14.6971 6.1882 -7.7474
17 6 2 0.5564 14.4784 5.3902 -7.2900
18 6 2 0.5564 15.2258 6.8450 -7.3017
19 7 1 -1.1128 -2.6192 10.3928 -0.0001
20 7 2 0.5564 -1.8647 10.9848 -0.5571
21 7 2 0.5564 -2.7265 9.2890 -0.0662
22 8 1 -1.1128 6.6045 2.0284 0.6044
23 8 2 0.5564 6.3655 1.0240 0.4799
24 8 2 0.5564 5.7494 2.3953 0.7463
25 9 1 -1.1128 8.8742 8.0839 4.6354
26 9 2 0.5564 9.7892 8.3767 4.4027
27 9 2 0.5564 8.3460 8.0813 3.8293
28 10 1 -1.1128 5.8097 5.6849 -0.7085
29 10 2 0.5564 5.9595 5.0452 -1.5572
30 10 2 0.5564 5.2079 6.4617 -0.9303
31 11 1 -1.1128 7.6988 2.7233 2.9973
32 11 2 0.5564 6.8675 3.0722 3.5888
33 11 2 0.5564 7.3299 2.4174 1.9170
34 12 1 -1.1128 11.0654 4.5876 3.5549
35 12 2 0.5564 10.3773 4.6342 2.8632
36 12 2 0.5564 10.4763 5.2035 4.2540
37 13 1 -1.1128 2.2356 2.9662 -0.3426
38 13 2 0.5564 2.8229 2.0247 -0.6116
39 13 2 0.5564 2.9868 3.0933 0.5014
40 14 1 -1.1128 0.6871 7.8853 -4.4203
41 14 2 0.5564 0.9105 7.8501 -3.3605
42 14 2 0.5564 0.0012 8.5407 -4.3902
43 15 1 -1.1128 6.1128 -0.4547 -7.5510
44 15 2 0.5564 6.0424 -0.9370 -6.6619
45 15 2 0.5564 6.9367 -0.2264 -7.6923
46 16 1 -1.1128 0.0269 7.9881 0.2323
47 16 2 0.5564 0.1935 7.5172 1.0657
48 16 2 0.5564 -1.0149 7.8781 0.1064
49 17 1 -1.1128 9.0979 4.4375 1.6117
50 17 2 0.5564 8.4785 3.5461 2.0370
51 17 2 0.5564 8.4603 5.1558 1.1598
52 18 1 -1.1128 8.0205 9.5154 -5.6604
53 18 2 0.5564 7.0891 9.9739 -5.8289
54 18 2 0.5564 8.2865 9.0556 -6.5545
55 19 1 -1.1128 7.1473 7.4240 -4.1023
56 19 2 0.5564 6.3979 6.9506 -4.4627
57 19 2 0.5564 7.1474 8.2756 -4.5903
58 20 1 -1.1128 4.2167 1.7495 -7.8409
59 20 2 0.5564 5.0855 1.2633 -7.6234
60 20 2 0.5564 3.7241 1.1296 -8.5781
61 21 1 -1.1128 5.4446 3.8853 -8.9068
62 21 2 0.5564 5.1557 4.7500 -8.6367
63 21 2 0.5564 4.8182 3.2069 -8.3794
64 22 1 -1.1128 6.0386 9.3204 -2.4511
65 22 2 0.5564 6.4152 8.5164 -2.8341
66 22 2 0.5564 6.7845 9.8389 -1.9444
67 23 1 -1.1128 0.1381 2.6793 4.5680
68 23 2 0.5564 -0.2809 3.4724 3.9858
69 23 2 0.5564 0.2112 2.0940 3.8931
70 24 1 -1.1128 11.3965 -3.5321 3.9345
71 24 2 0.5564 11.7433 -3.3427 4.8119
72 24 2 0.5564 11.9756 -4.4481 3.3455
73 25 1 -1.1128 5.1961 -2.3263 -5.0160
74 25 2 0.5564 5.5845 -2.6326 -3.9068
75 25 2 0.5564 4.9080 -1.3876 -5.0752
76 26 1 -1.1128 1.9884 4.6024 -2.2437
77 26 2 0.5564 1.1159 4.2702 -2.6316
78 26 2 0.5564 2.1782 3.8946 -1.7102
79 27 1 -1.1128 3.4286 -3.8254 2.9948
80 27 2 0.5564 4.1401 -3.1922 2.7029
81 27 2 0.5564 2.6837 -3.3946 2.5442
82 28 1 -1.1128 3.4972 4.0312 -4.2534
83 28 2 0.5564 3.2969 3.6780 -5.1881
84 28 2 0.5564 2.8283 3.9277 -3.5286
85 29 1 -1.1128 2.5581 3.6391 -6.8719
86 29 2 0.5564 3.0211 2.9680 -7.2340
87 29 2 0.5564 1.6628 3.1519 -7.0428
88 30 1 -1.1128 5.1229 6.3022 -5.5820
89 30 2 0.5564 4.4630 7.1815 -5.7854
90 30 2 0.5564 4.3184 5.7177 -5.1800
91 31 1 -1.1128 7.5284 -1.4968 -1.3872
92 31 2 0.5564 7.6646 -0.7307 -2.1238
93 31 2 0.5564 8.3855 -1.4659 -1.0203
94 32 1 -1.1128 8.0180 3.6712 -5.5542
95 32 2 0.5564 8.2173 4.4555 -6.0837
96 32 2 0.5564 8.6346 3.0360 -5.6917
97 33 1 -1.1128 16.0706 -0.2590 -5.0003
98 33 2 0.5564 16.0528 -0.2534 -5.8842
99 33 2 0.5564 15.1342 -0.1953 -4.7488
100 34 1 -1.1128 8.7867 0.2163 4.2394
101 34 2 0.5564 8.4837 1.0587 3.5818
102 34 2 0.5564 9.5314 -0.1659 3.6191
103 35 1 -1.1128 3.2350 8.1868 -5.7667
104 35 2 0.5564 2.4833 8.4024 -5.2868
105 35 2 0.5564 3.9304 8.8555 -5.5968
106 36 1 -1.1128 9.4233 5.7040 5.6841
107 36 2 0.5564 10.0973 5.7599 6.4736
108 36 2 0.5564 8.9638 6.5622 5.3889
109 37 1 -1.1128 5.4246 -2.0137 2.6572
110 37 2 0.5564 5.8998 -1.6969 3.5099
111 37 2 0.5564 5.5652 -1.5574 1.8565
112 38 1 -1.1128 -0.1317 1.9336 0.3190
113 38 2 0.5564 0.7894 2.2269 -0.1673
114 38 2 0.5564 -0.9403 2.4988 0.0245
115 39 1 -1.1128 11.8227 3.8276 -2.8190
116 39 2 0.5564 11.9045 2.9412 -3.1683
117 39 2 0.5564 11.4467 3.8091 -2.0650
118 40 1 -1.1128 -1.2368 9.7290 -3.9382
119 40 2 0.5564 -1.9257 9.0012 -3.6423
120 40 2 0.5564 -1.6406 10.5300 -4.3785
121 41 1 -1.1128 10.0754 7.7161 -0.3734
122 41 2 0.5564 10.0574 7.7060 -1.5087
123 41 2 0.5564 9.2537 7.2615 -0.1815
124 42 1 -1.1128 -1.0123 6.0538 -4.7618
125 42 2 0.5564 -0.8778 5.3177 -4.1260
126 42 2 0.5564 -0.2149 6.6635 -4.6539
127 43 1 -1.1128 -1.7545 11.2259 2.6532
128 43 2 0.5564 -1.2454 10.3615 3.0272
129 43 2 0.5564 -2.3173 10.8721 1.8706
130 44 1 -1.1128 12.4899 11.0304 5.7156
131 44 2 0.5564 11.9240 11.7099 5.4868
132 44 2 0.5564 13.3664 11.1913 5.4911
133 45 1 -1.1128 9.6108 7.3259 -3.0932
134 45 2 0.5564 8.6400 7.3761 -3.4426
135 45 2 0.5564 10.1122 6.7757 -3.8567
136 46 1 -1.1128 2.6948 11.0492 4.8330
137 46 2 0.5564 2.7024 11.5790 4.0320
138 46 2 0.5564 3.0263 9.9915 4.5693
139 47 1 -1.1128 2.3981 10.4891 -1.9366
140 47 2 0.5564 3.2686 9.9508 -1.8156
141 47 2 0.5564 2.0376 10.5907 -0.9310
142 48 1 -1.1128 0.9756 10.6255 0.5327
143 48 2 0.5564 0.1537 11.3288 0.6051
144 48 2 0.5564 0.5959 9.7905 0.6015
145 49 1 -1.1128 6.8518 4.4218 -2.9838
146 49 2 0.5564 7.5360 4.5692 -4.0439
147 49 2 0.5564 6.2902 3.8048 -3.2906
148 50 1 -1.1128 5.2164 6.6410 4.0581
149 50 2 0.5564 5.2488 6.3928 4.9660
150 50 2 0.5564 4.2955 7.0350 3.8882
151 51 1 -1.1128 0.3549 6.8498 2.9794
152 51 2 0.5564 1.0995 6.6230 3.4814
153 51 2 0.5564 -0.3406 6.0721 3.1076
154 52 1 -1.1128 8.2297 2.5219 -1.4475
155 52 2 0.5564 7.9524 3.3343 -1.9877
156 52 2 0.5564 7.7019 2.4769 -0.5961
157 53 1 -1.1128 7.8245 6.7224 0.7167
158 53 2 0.5564 7.2205 6.3967 0.2032
159 53 2 0.5564 7.5432 7.1690 1.5959
160 54 1 -1.1128 10.4529 3.5595 -0.4618
161 54 2 0.5564 9.6554 3.0113 -1.0854
162 54 2 0.5564 9.9860 3.8155 0.4575
163 55 1 -1.1128 1.6091 7.3889 -1.9372
164 55 2 0.5564 0.8772 7.6442 -1.3852
165 55 2 0.5564 1.6148 6.3199 -2.0961
166 56 1 -1.1128 4.5847 3.7052 0.8854
167 56 2 0.5564 4.9276 4.3327 0.4075
168 56 2 0.5564 4.7321 3.6932 1.9221
169 57 1 -1.1128 7.7489 0.5850 -3.1580
170 57 2 0.5564 8.1265 1.3116 -2.5648
171 57 2 0.5564 6.8161 0.9453 -3.2322
172 58 1 -1.1128 1.0484 11.8440 -4.0472
173 58 2 0.5564 0.7280 11.1323 -4.6797
174 58 2 0.5564 1.1330 11.2806 -3.2996
175 59 1 -1.1128 4.8434 1.5791 -3.4338
176 59 2 0.5564 4.6193 2.5111 -3.8020
177 59 2 0.5564 4.5590 0.9918 -4.2585
178 60 1 -1.1128 9.5745 -0.8244 7.3550
179 60 2 0.5564 9.1589 -1.6054 6.9407
180 60 2 0.5564 8.9480 -0.4525 7.9368
181 61 1 -1.1128 10.2730 1.8899 6.5393
182 61 2 0.5564 10.8570 2.0184 5.8089
183 61 2 0.5564 9.9510 0.9418 6.2056
184 62 1 -1.1128 4.0753 7.8192 -0.7739
185 62 2 0.5564 4.6355 8.0586 -1.5322
186 62 2 0.5564 3.2914 7.7124 -1.2255
187 63 1 -1.1128 3.7702 0.4562 -1.1286
188 63 2 0.5564 3.4391 -0.4159 -1.3938
189 63 2 0.5564 4.1240 0.8020 -1.9818
190 64 1 -1.1128 6.9868 8.1387 2.7851
191 64 2 0.5564 6.2262 7.6175 3.5031
192 64 2 0.5564 6.8162 9.0302 2.5521
Bonds
1 1 1 2
2 1 1 3
3 1 4 5
4 1 4 6
5 1 7 8
6 1 7 9
7 1 10 11
8 1 10 12
9 1 13 14
10 1 13 15
11 1 16 17
12 1 16 18
13 1 19 20
14 1 19 21
15 1 22 23
16 1 22 24
17 1 25 26
18 1 25 27
19 1 28 29
20 1 28 30
21 1 31 32
22 1 31 33
23 1 34 35
24 1 34 36
25 1 37 38
26 1 37 39
27 1 40 41
28 1 40 42
29 1 43 44
30 1 43 45
31 1 46 47
32 1 46 48
33 1 49 50
34 1 49 51
35 1 52 53
36 1 52 54
37 1 55 56
38 1 55 57
39 1 58 59
40 1 58 60
41 1 61 62
42 1 61 63
43 1 64 65
44 1 64 66
45 1 67 68
46 1 67 69
47 1 70 71
48 1 70 72
49 1 73 74
50 1 73 75
51 1 76 77
52 1 76 78
53 1 79 80
54 1 79 81
55 1 82 83
56 1 82 84
57 1 85 86
58 1 85 87
59 1 88 89
60 1 88 90
61 1 91 92
62 1 91 93
63 1 94 95
64 1 94 96
65 1 97 98
66 1 97 99
67 1 100 101
68 1 100 102
69 1 103 104
70 1 103 105
71 1 106 107
72 1 106 108
73 1 109 110
74 1 109 111
75 1 112 113
76 1 112 114
77 1 115 116
78 1 115 117
79 1 118 119
80 1 118 120
81 1 121 122
82 1 121 123
83 1 124 125
84 1 124 126
85 1 127 128
86 1 127 129
87 1 130 131
88 1 130 132
89 1 133 134
90 1 133 135
91 1 136 137
92 1 136 138
93 1 139 140
94 1 139 141
95 1 142 143
96 1 142 144
97 1 145 146
98 1 145 147
99 1 148 149
100 1 148 150
101 1 151 152
102 1 151 153
103 1 154 155
104 1 154 156
105 1 157 158
106 1 157 159
107 1 160 161
108 1 160 162
109 1 163 164
110 1 163 165
111 1 166 167
112 1 166 168
113 1 169 170
114 1 169 171
115 1 172 173
116 1 172 174
117 1 175 176
118 1 175 177
119 1 178 179
120 1 178 180
121 1 181 182
122 1 181 183
123 1 184 185
124 1 184 186
125 1 187 188
126 1 187 189
127 1 190 191
128 1 190 192
Angles
1 1 2 1 3
2 1 5 4 6
3 1 8 7 9
4 1 11 10 12
5 1 14 13 15
6 1 17 16 18
7 1 20 19 21
8 1 23 22 24
9 1 26 25 27
10 1 29 28 30
11 1 32 31 33
12 1 35 34 36
13 1 38 37 39
14 1 41 40 42
15 1 44 43 45
16 1 47 46 48
17 1 50 49 51
18 1 53 52 54
19 1 56 55 57
20 1 59 58 60
21 1 62 61 63
22 1 65 64 66
23 1 68 67 69
24 1 71 70 72
25 1 74 73 75
26 1 77 76 78
27 1 80 79 81
28 1 83 82 84
29 1 86 85 87
30 1 89 88 90
31 1 92 91 93
32 1 95 94 96
33 1 98 97 99
34 1 101 100 102
35 1 104 103 105
36 1 107 106 108
37 1 110 109 111
38 1 113 112 114
39 1 116 115 117
40 1 119 118 120
41 1 122 121 123
42 1 125 124 126
43 1 128 127 129
44 1 131 130 132
45 1 134 133 135
46 1 137 136 138
47 1 140 139 141
48 1 143 142 144
49 1 146 145 147
50 1 149 148 150
51 1 152 151 153
52 1 155 154 156
53 1 158 157 159
54 1 161 160 162
55 1 164 163 165
56 1 167 166 168
57 1 170 169 171
58 1 173 172 174
59 1 176 175 177
60 1 179 178 180
61 1 182 181 183
62 1 185 184 186
63 1 188 187 189
64 1 191 190 192

View File

@ -0,0 +1,29 @@
units real
atom_style gle # 6
pair_style lj/cut/tip4p/long 1 2 1 1 0.14714951 8
bond_style class2
angle_style harmonic
kspace_style pppm 0.0001
read_data data.tip4pf
pair_coeff * 2 0 0
pair_coeff 1 1 0.1852 3.1589022
#q-tip4pf bond parameters
bond_coeff 1 0.9419 607.19354 -1388.6516 1852.577
angle_coeff 1 43.93 107.4
thermo 100
thermo_style custom step temp pe ke etotal
timestep 0.5
dump snapxyz all atom 100 h2o.lammpstraj
# some problem
fix 1 all gle 6 smart.A 300 31415
fix_modify 1 energy yes
#fix 1 all nve
run 10000

View File

@ -19,7 +19,7 @@ PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
user-cuda user-eff user-fep user-lb user-misc user-molfile \ user-cuda user-eff user-fep user-lb user-misc user-molfile \
user-omp user-phonon user-qmmm user-reaxc user-sph user-omp user-phonon user-qmmm user-reaxc user-sph user-gle
PACKLIB = gpu kim meam poems reax voronoi \ PACKLIB = gpu kim meam poems reax voronoi \
user-atc user-awpmd user-colvars user-qmmm user-cuda user-molfile user-atc user-awpmd user-colvars user-qmmm user-cuda user-molfile

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,99 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
* A class for an atom kind which carries around a variable number of
* additional degrees of freedom, to integrate a Generalized Langevin
* Equation dynamics.
*
* TODO: Probably this is better implemented as a hybrid type.
* ---------------------------------------------------------------------- */
#ifdef ATOM_CLASS
AtomStyle(gle,AtomVecGLE)
#else
#ifndef LMP_ATOM_VEC_GLE_H
#define LMP_ATOM_VEC_GLE_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecGLE : public AtomVec {
public:
AtomVecGLE(class LAMMPS *); //, int, char **); AtomVec classes used to accept arguments
virtual ~AtomVecGLE() {}
void grow(int);
void grow_reset();
void copy(int, int, int);
virtual int pack_comm(int, int *, double *, int, int *);
virtual int pack_comm_vel(int, int *, double *, int, int *);
virtual void unpack_comm(int, int, double *);
virtual void unpack_comm_vel(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
virtual int pack_border(int, int *, double *, int, int *);
virtual int pack_border_vel(int, int *, double *, int, int *);
int pack_border_hybrid(int, int *, double *);
virtual void unpack_border(int, int, double *);
virtual void unpack_border_vel(int, int, double *);
int unpack_border_hybrid(int, int, double *);
virtual int pack_exchange(int, double *);
virtual int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, imageint, char **);
int data_atom_hybrid(int, char **);
void pack_data(double **);
int pack_data_hybrid(int, double *);
void write_data(FILE *, int, double **);
int write_data_hybrid(FILE *, double *);
bigint memory_usage();
protected:
int ns;
tagint *tag;
int *type,*mask;
imageint *image;
// s[i][j] i is the atom index, j contains 3*(ns+1) items
// corresponding to x,y,z coordinates of different
// additional DOFs, i.e. s[0] will be [ x0 y0 z0 x1 y1 z1 ... ]
double **x, **v, **f, **s;
double *q;
tagint *molecule;
int **nspecial;
tagint **special;
int *num_bond;
int **bond_type;
tagint **bond_atom;
int *num_angle;
int **angle_type;
tagint **angle_atom1,**angle_atom2,**angle_atom3;
int *num_dihedral;
int **dihedral_type;
tagint **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4;
int *num_improper;
int **improper_type;
tagint **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
};
}
#endif
#endif

522
src/USER-GLE/fix_gle.cpp Normal file
View File

@ -0,0 +1,522 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Michele Ceriotti (Oxford), Joe Morrone (Columbia)
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_gle.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "comm.h"
#include "input.h"
#include "variable.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "group.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{NOBIAS,BIAS};
enum{CONSTANT,EQUAL,ATOM};
//#define GLE_DEBUG 1
/* syntax for fix_gle:
* fix_gle nfix id-group gle ns amatrix [temp|cmatrix] seed
*
* */
/* ---------------------------------------------------------------------- */
namespace GLE {
/*!!MC!! here begins the GLE stuff*/
#define midx(n,i,j) ((i)*(n)+(j))
//"stabilized" cholesky decomposition. does a LDL^t decomposition, then sets to zero the negative diagonal elements and gets MM^t
void StabCholesky(unsigned long n, const double* MMt, double* M)
{
double L[n*n], D[n];
unsigned long i,j,k;
for(i=0; i<n; ++i) D[i]=0.;
for(i=0; i<n*n; ++i) L[i]=0.;
for(i=0; i<n; ++i)
{
L[midx(n,i,i)]=1.;
for (j=0; j<i; j++)
{
printf("%d %d\n", i,j);
L[midx(n,i,j)]=MMt[midx(n,i,j)];
for (k=0; k<j; ++k) L[midx(n,i,j)]-=L[midx(n,i,k)]*L[midx(n,j,k)]*D[k];
if (D[j]!=0.) L[midx(n,i,j)]/=D[j];
else L[midx(n,i,j)]=0.0;
}
D[i]=MMt[midx(n,i,i)];
for (k=0; k<i; ++k) D[i]-=L[midx(n,i,k)]*L[midx(n,i,k)]*D[k];
}
for(i=0; i<n; ++i) D[i]=(D[i]>0.?sqrt(D[i]):0.);
for(i=0; i<n; ++i) for (j=0; j<n; j++) M[midx(n,i,j)]=L[midx(n,i,j)]*D[j];
}
void MyMult(unsigned long n, unsigned long m, unsigned long r, const double* A, const double* B, double* C, double cf=0.0)
{
// !! TODO !! should probably call BLAS (or check if some efficient matrix-matrix multiply is implemented somewhere in LAMMPS)
// (rows x cols) :: A is n x r, B is r x m, C is n x m
unsigned long i,j,k; double *cij;
for (i=0; i<n; ++i)
for (j=0; j<m; ++j)
{
cij=&C[midx(m,i,j)]; *cij *= cf;
for (k=0; k<r; ++k) *cij+=A[midx(r,i,k)]*B[midx(m,k,j)];
}
}
void MyTrans(unsigned long n, const double* A, double* AT)
{
for (unsigned long i=0; i<n; ++i) for (unsigned long j=0; j<n; ++j) AT[j*n+i]=A[i*n+j];
}
void MyPrint(unsigned long n, const double* A)
{
for (unsigned long k=0; k<n*n; ++k) { printf("%10.5e ", A[k]); if ((k+1)%n==0) printf("\n");}
}
//matrix exponential by scaling and squaring.
void MatrixExp(unsigned long n, const double* M, double* EM, unsigned long j=8, unsigned long k=8)
{
double tc[j+1], SM[n*n], TMP[n*n];
double onetotwok=pow(0.5,1.0*k);
tc[0]=1;
for (int i=0; i<j; ++i) tc[i+1]=tc[i]/(i+1.0);
for (int i=0; i<n*n; ++i) { SM[i]=M[i]*onetotwok; EM[i]=0.0; TMP[i]=0.0; }
for (int i=0; i<n; ++i) EM[midx(n,i,i)]=tc[j];
//taylor exp of scaled matrix
for (int p=j-1; p>=0; p--)
{
MyMult(n, n, n, SM, EM, TMP); for (int i=0; i<n*n; ++i) EM[i]=TMP[i];
for (int i=0; i<n; ++i) EM[midx(n,i,i)]+=tc[p];
}
for (int p=0; p<k; p++)
{ MyMult(n, n, n, EM, EM, TMP); for (int i=0; i<n*n; ++i) EM[i]=TMP[i]; }
}
}
FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 5) error->all(FLERR,"Illegal fix gle command");
ns = atof(arg[3]);
if (ns>atom->ns) error->all(FLERR, "Atom kind has not enough GLE slots");
A = new double[(ns+1)*(ns+1)];
C = new double[(ns+1)*(ns+1)];
T=S=gle_tmp=gle_rnd=NULL; gle_buff=0;
// LOADS A matrix
printf("Reading A from %s\n", arg[4]);
FILE* fgle = fopen(arg[4], "r");
if (fgle==NULL) error->all(FLERR, "Cannot open A matrix GLE input");
for (int i=0; i<(ns+1)*(ns+1); ++i) //ASSUMES A IS IN THE CORRECT INVERSE TIME UNITS
if (!fscanf(fgle,"%lf",&(A[i]))) error->all(FLERR, "Cannot read A matrix GLE input");
fclose(fgle);
temp=atof(arg[5]);
if (temp==0.0)
{
printf("Reading C from %s\n", arg[5]);
fgle = fopen(arg[5], "r");
if (fgle==NULL) error->all(FLERR, "Cannot open C matrix GLE input");
for (int i=0; i<(ns+1)*(ns+1); ++i) //ASSUMES C IS IN THE CORRECT TEMPERATURE UNITS
if (!fscanf(fgle,"%lf",&(C[i]))) error->all(FLERR, "Cannot read C matrix GLE input");
}
else
{
for (int i=0; i<(ns+1)*(ns+1); ++i) C[i]=0.0;
for (int i=0; i<(ns+1)*(ns+1); i+=(ns+2)) C[i]=temp*force->boltz/force->mvv2e;
}
#ifdef GLE_DEBUG
printf("A Matrix\n");
GLE::MyPrint(ns+1,A);
printf("C Matrix\n");
GLE::MyPrint(ns+1,C);
#endif
tau=10;
int seed = atoi(arg[6]);
// initialize Marsaglia RNG with processor-unique seed
// NB: this means runs will not be the same with different numbers of processors
random = new RanMars(lmp,seed + comm->me);
if (seed <= 0) error->all(FLERR,"Illegal fix gle command");
// allocate per-type arrays for mass-scaling
sqrt_m = new double[atom->ntypes+1];
energy = 0.0;
}
/* ---------------------------------------------------------------------- */
FixGLE::~FixGLE()
{
delete random;
delete [] sqrt_m;
delete [] A; delete [] C; delete [] S; delete [] T;
delete [] gle_tmp; delete [] gle_rnd;
}
/* ---------------------------------------------------------------------- */
int FixGLE::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
mask |= THERMO_ENERGY;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixGLE::init()
{
dogle = 1;
dorattle = 1;
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
// compute Langevin terms
T = new double[(ns+1)*(ns+1)]; S = new double[(ns+1)*(ns+1)];
double *tmp1 = new double[(ns+1)*(ns+1)], *tmp2 = new double[(ns+1)*(ns+1)];
for (int i=0; i<(ns+1)*(ns+1); ++i) { tmp1[i]=-A[i]*update->dt*0.5; tmp2[i]=S[i]=0.0; }
GLE::MatrixExp(ns+1,tmp1,T);
GLE::MyMult(ns+1,ns+1,ns+1,T,C,tmp1);
GLE::MyTrans(ns+1,T,tmp2);
GLE::MyMult(ns+1,ns+1,ns+1,tmp1,tmp2,S);
for (int i=0; i<(ns+1)*(ns+1); ++i) tmp1[i]=C[i]-S[i];
GLE::StabCholesky(ns+1, tmp1, S); //!TODO use symmetric square root, which is more stable.
#ifdef GLE_DEBUG
printf("T Matrix\n");
GLE::MyPrint(ns+1,T);
printf("S Matrix\n");
GLE::MyPrint(ns+1,S);
#endif
delete[] tmp1; delete[] tmp2;
// set force prefactors
if (!atom->rmass) {
for (int i = 1; i <= atom->ntypes; i++) {
sqrt_m[i] = sqrt(atom->mass[i]);
}
}
if (strstr(update->integrate_style,"respa")) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
}
}
/* ---------------------------------------------------------------------- */
void FixGLE::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
void FixGLE::gle_integrate()
{
double **v = atom->v;
double **s = atom->s;
double *rmass = atom->rmass, smi, ismi;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (nlocal > gle_buff)
{
if (gle_tmp !=NULL) {
delete [] gle_tmp; delete [] gle_rnd;
}
gle_tmp = new double[nlocal*(ns+1)];
gle_rnd = new double[nlocal*(ns+1)];
gle_buff = nlocal;
}
#ifdef GLE_DEBUG
printf("!MC! GLE THERMO STEP dt=%f\n",update->dt);
#endif
// gets kinetic energy before doing anything to the velocities
double deltae=0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ismi=rmass[i];
for (int j = 0; j<3; ++j) { deltae-=ismi*v[i][j]*v[i][j]; }
} }
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ismi=atom->mass[type[i]];
for (int j = 0; j<3; ++j) { deltae-=ismi*v[i][j]*v[i][j]; }
} }
}
// s is just taken to be the mass-scaled velocity!
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
smi=sqrt(rmass[i]);
for (int j = 0; j<3; ++j) { s[i][j] = v[i][j]*smi; }
} }
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
smi=sqrt_m[type[i]];
for (int j = 0; j<3; ++j) { s[i][j] = v[i][j]*smi; }
} }
}
for (int k = 0; k<3; k++) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) for (int j = 0; j<(ns+1); ++j) gle_rnd[j*(nlocal)+i] = s[i][3*j+k];
}
GLE::MyMult(ns+1,nlocal,ns+1,T,gle_rnd,gle_tmp);
for (int i = 0; i < nlocal*(ns+1); i++) gle_rnd[i] = random->gaussian();
GLE::MyMult(ns+1,nlocal,ns+1,S,gle_rnd,gle_tmp,1.0);
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) for (int j = 0; j<(ns+1); ++j) s[i][3*j+k] = gle_tmp[j*(nlocal)+i];
}
}
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ismi=1.0/sqrt(rmass[i]);
for (int j = 0; j<3; ++j) { v[i][j] = s[i][j]*ismi; }
} }
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ismi=1.0/sqrt_m[type[i]];
for (int j = 0; j<3; ++j) { v[i][j] = s[i][j]*ismi; }
} }
}
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ismi=rmass[i];
for (int j = 0; j<3; ++j) { deltae+=ismi*v[i][j]*v[i][j]; }
} }
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ismi=atom->mass[type[i]];
for (int j = 0; j<3; ++j) { deltae+=ismi*v[i][j]*v[i][j]; }
} }
}
energy += deltae*0.5*force->mvv2e;
}
void FixGLE::initial_integrate(int vflag)
{
double dtfm;
// update v and x of atoms in group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (dogle) gle_integrate();
#ifdef GLE_DEBUG
printf("!MC! GLE P1 STEP dt=%f\n",dtv);
printf("!MC! GLE Q STEP\n");
#endif
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
void FixGLE::final_integrate()
{
double dtfm;
// update v of atoms in group
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
#ifdef GLE_DEBUG
printf("!MC! GLE P2 STEP dt=%f\n",dtv);
#endif
if (dogle) gle_integrate();
}
/* ---------------------------------------------------------------------- */
void FixGLE::initial_integrate_respa(int vflag, int ilevel, int iloop)
{
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
// innermost level - NVE update of v and x
// all other levels - NVE update of v
if (ilevel==nlevels_respa-1) gle_integrate();
dogle=0;
if (ilevel == 0) initial_integrate(vflag);
else { dorattle=1; final_integrate();} // If RATTLE S2 should not be called here, dorattle=0
}
void FixGLE::final_integrate_respa(int ilevel, int iloop)
{
dtv = step_respa[ilevel]; //!MC!
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dogle=0;
if(0) { // If true turn on stage 2 rattle for final integrate only
if(ilevel==nlevels_respa-1)
dorattle=1;
else
dorattle=0;
} else { dorattle=1; }
final_integrate();
if (ilevel==nlevels_respa-1) gle_integrate();
}
double FixGLE::compute_scalar()
{
double energy_me = energy;
double energy_all;
MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
return -energy_all;
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */
void *FixGLE::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"t_target") == 0) {
return &temp;
}
return NULL;
}

61
src/USER-GLE/fix_gle.h Normal file
View File

@ -0,0 +1,61 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(gle,FixGLE)
#else
#ifndef LMP_FIX_GLE_H
#define LMP_FIX_GLE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixGLE : public Fix {
public:
FixGLE(class LAMMPS *, int, char **);
virtual ~FixGLE();
int setmask();
void init();
void setup(int);
void gle_integrate();
void initial_integrate_respa(int vflag, int ilevel, int iloop);
void final_integrate_respa(int ilevel, int iloop);
void initial_integrate(int vflag);
void final_integrate();
double compute_scalar();
virtual void *extract(const char *, int &);
protected:
int ns;
double *A, *C, *S, *T;
double *gle_tmp, *gle_rnd; int gle_buff;
double temp, tau, dtv, dtf;
int dogle;
int dorattle;
class RanMars *random;
double *sqrt_m;
double *step_respa;
double energy;
int nlevels_respa;
double **vaux;
};
}
#endif
#endif

View File

@ -52,6 +52,7 @@ using namespace MathConst;
Atom::Atom(LAMMPS *lmp) : Pointers(lmp) Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
{ {
natoms = 0; natoms = 0;
ns = 0; //!GLE
nlocal = nghost = nmax = 0; nlocal = nghost = nmax = 0;
ntypes = 0; ntypes = 0;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
@ -72,6 +73,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
type = mask = NULL; type = mask = NULL;
image = NULL; image = NULL;
x = v = f = NULL; x = v = f = NULL;
s = NULL; //!GLE
molecule = NULL; molecule = NULL;
molindex = molatom = NULL; molindex = molatom = NULL;
@ -205,6 +207,7 @@ Atom::~Atom()
memory->destroy(x); memory->destroy(x);
memory->destroy(v); memory->destroy(v);
memory->destroy(f); memory->destroy(f);
memory->destroy(s); //!GLE
memory->destroy(molecule); memory->destroy(molecule);
memory->destroy(molindex); memory->destroy(molindex);
@ -1891,6 +1894,7 @@ void *Atom::extract(char *name)
if (strcmp(name,"x") == 0) return (void *) x; if (strcmp(name,"x") == 0) return (void *) x;
if (strcmp(name,"v") == 0) return (void *) v; if (strcmp(name,"v") == 0) return (void *) v;
if (strcmp(name,"f") == 0) return (void *) f; if (strcmp(name,"f") == 0) return (void *) f;
if (strcmp(name,"s") == 0) return (void *) s; //!GLE
if (strcmp(name,"molecule") == 0) return (void *) molecule; if (strcmp(name,"molecule") == 0) return (void *) molecule;
if (strcmp(name,"q") == 0) return (void *) q; if (strcmp(name,"q") == 0) return (void *) q;
if (strcmp(name,"mu") == 0) return (void *) mu; if (strcmp(name,"mu") == 0) return (void *) mu;

View File

@ -51,6 +51,7 @@ class Atom : protected Pointers {
imageint *image; imageint *image;
double **x,**v,**f; double **x,**v,**f;
double **s; int ns; //!GLE
tagint *molecule; tagint *molecule;
int *molindex,*molatom; int *molindex,*molatom;