Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionAtomistic.h"
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/Vector.h"
26 : #include "tools/Matrix.h"
27 : #include "tools/AtomNumber.h"
28 : #include "tools/Tools.h"
29 : #include "core/Atoms.h"
30 : #include "tools/Pbc.h"
31 :
32 : #include <vector>
33 : #include <string>
34 :
35 : using namespace std;
36 :
37 : namespace PLMD {
38 : namespace generic {
39 :
40 : //+PLUMEDOC GENERIC RESET_CELL
41 : /*
42 : This action is used to rotate the full cell
43 :
44 : This can be used to modify the periodic box. Notice that
45 : this is done at fixed scaled coordinates,
46 : so that also atomic coordinates for the entire system are affected.
47 : To see what effect try
48 : the \ref DUMPATOMS directive to output the atomic positions.
49 :
50 : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
51 : after rotation. See also \ref FIT_TO_TEMPLATE
52 :
53 : Currently, only TYPE=TRIANGULAR is implemented, which allows one to reset
54 : the cell to a lower triangular one. Namely, a proper rotation is found that allows
55 : rotating the box so that the first lattice vector is in the form (ax,0,0),
56 : the second lattice vector is in the form (bx,by,0), and the third lattice vector is
57 : arbitrary.
58 :
59 : \attention
60 : The implementation of this action is available but should be considered in testing phase. Please report any
61 : strange behavior.
62 :
63 : \attention
64 : This directive modifies the stored position at the precise moment
65 : it is executed. This means that only collective variables
66 : which are below it in the input script will see the corrected positions.
67 : Unless you
68 : know exactly what you are doing, leave the default stride (1), so that
69 : this action is performed at every MD step.
70 :
71 : \par Examples
72 :
73 : Reset cell to be triangular after a rototranslational fit
74 : \verbatim
75 : DUMPATOMS FILE=dump-original.xyz ATOMS=1-20
76 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=OPTIMAL
77 : DUMPATOMS FILE=dump-fit.xyz ATOMS=1-20
78 : RESET_CELL TYPE=TRIANGULAR
79 : DUMPATOMS FILE=dump-reset.xyz ATOMS=1-20
80 : \endverbatim
81 : (see also \ref DUMPATOMS)
82 :
83 :
84 : */
85 : //+ENDPLUMEDOC
86 :
87 :
88 4 : class ResetCell:
89 : public ActionPilot,
90 : public ActionAtomistic
91 : {
92 : std::string type;
93 : Tensor rotation,newbox;
94 :
95 : public:
96 : explicit ResetCell(const ActionOptions&ao);
97 : static void registerKeywords( Keywords& keys );
98 : void calculate();
99 : void apply();
100 : };
101 :
102 2525 : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
103 :
104 3 : void ResetCell::registerKeywords( Keywords& keys ) {
105 3 : Action::registerKeywords( keys );
106 3 : ActionAtomistic::registerKeywords( keys );
107 3 : keys.add("compulsory","STRIDE","1","the frequency with which molecules are reassembled. Unless you are completely certain about what you are doing leave this set equal to 1!");
108 3 : keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
109 3 : }
110 :
111 2 : ResetCell::ResetCell(const ActionOptions&ao):
112 : Action(ao),
113 : ActionPilot(ao),
114 2 : ActionAtomistic(ao)
115 : {
116 2 : type.assign("TRIANGULAR");
117 2 : parse("TYPE",type);
118 :
119 2 : log<<" type: "<<type<<"\n";
120 2 : if(type=="TRIANGULAR") {
121 0 : } else error("undefined type "+type);
122 :
123 2 : checkRead();
124 2 : }
125 :
126 :
127 17 : void ResetCell::calculate() {
128 :
129 17 : Pbc & pbc(modifyGlobalPbc());
130 :
131 17 : Tensor box=pbc.getBox();
132 :
133 : // moduli of lattice vectors
134 17 : double a=modulo(box.getRow(0));
135 17 : double b=modulo(box.getRow(1));
136 17 : double c=modulo(box.getRow(2));
137 : // cos-angle between lattice vectors
138 17 : double ab=dotProduct(box.getRow(0),box.getRow(1))/(a*b);
139 17 : double ac=dotProduct(box.getRow(0),box.getRow(2))/(a*c);
140 17 : double bc=dotProduct(box.getRow(1),box.getRow(2))/(b*c);
141 :
142 : // generate a new set of lattice vectors as a lower triangular matrix
143 17 : newbox[0][0]=a;
144 17 : newbox[1][0]=b*ab;
145 17 : newbox[1][1]=std::sqrt(b*b-newbox[1][0]*newbox[1][0]);
146 17 : newbox[2][0]=c*ac;
147 17 : newbox[2][1]=c*(bc-ac*ab)/std::sqrt(1-ab*ab);
148 17 : newbox[2][2]=std::sqrt(c*c-newbox[2][0]*newbox[2][0]-newbox[2][1]*newbox[2][1]);
149 :
150 17 : if(determinant(newbox)*determinant(box)<0) newbox[2][2]=-newbox[2][2];
151 :
152 : // rotation matrix from old to new coordinates
153 17 : rotation=transpose(matmul(inverse(box),newbox));
154 :
155 : // rotate all coordinates
156 1611 : for(unsigned i=0; i<getTotAtoms(); i++) {
157 1594 : Vector & ato (modifyPosition(AtomNumber::index(i)));
158 1594 : ato=matmul(rotation,ato);
159 : }
160 : // rotate box
161 17 : pbc.setBox(newbox);
162 17 : }
163 :
164 17 : void ResetCell::apply() {
165 : // rotate back forces
166 1611 : for(unsigned i=0; i<getTotAtoms(); i++) {
167 1594 : Vector & f(modifyGlobalForce(AtomNumber::index(i)));
168 1594 : f=matmul(transpose(rotation),f);
169 : }
170 :
171 17 : Tensor& virial(modifyGlobalVirial());
172 : // I have no mathematical derivation for this.
173 : // The reasoning is the following.
174 : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
175 : // The final virial should be rotationally invariant, that is symmetric.
176 : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
177 : // be changed so as to enforce rotational invariance. Thus we here have to
178 : // make the virial matrix symmetric.
179 : // Since h^T is upper triangular, it can be shown that any change in these elements
180 : // will only affect the corresponding elements of the virial matrix.
181 : // Thus, the only possibility is to set the corresponding elements
182 : // of the virial matrix equal to their symmetric ones.
183 : // GB
184 17 : virial[0][1]=virial[1][0];
185 17 : virial[0][2]=virial[2][0];
186 17 : virial[1][2]=virial[2][1];
187 : // rotate back virial
188 17 : virial=matmul(transpose(rotation),matmul(virial,rotation));
189 :
190 :
191 :
192 17 : }
193 :
194 : }
195 2523 : }
|
