Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 : \plumedfile
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 : \endplumedfile
81 :
82 :
83 : */
84 : //+ENDPLUMEDOC
85 :
86 :
87 8 : class ResetCell:
88 : public ActionPilot,
89 : public ActionAtomistic
90 : {
91 : std::string type;
92 : Tensor rotation,newbox;
93 :
94 : public:
95 : explicit ResetCell(const ActionOptions&ao);
96 : static void registerKeywords( Keywords& keys );
97 : void calculate();
98 : void apply();
99 : };
100 :
101 6454 : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
102 :
103 3 : void ResetCell::registerKeywords( Keywords& keys ) {
104 3 : Action::registerKeywords( keys );
105 3 : ActionAtomistic::registerKeywords( keys );
106 15 : 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!");
107 15 : keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
108 3 : }
109 :
110 2 : ResetCell::ResetCell(const ActionOptions&ao):
111 : Action(ao),
112 : ActionPilot(ao),
113 4 : ActionAtomistic(ao)
114 : {
115 2 : type.assign("TRIANGULAR");
116 4 : parse("TYPE",type);
117 :
118 2 : log<<" type: "<<type<<"\n";
119 2 : if(type=="TRIANGULAR") {
120 0 : } else error("undefined type "+type);
121 :
122 2 : checkRead();
123 2 : }
124 :
125 :
126 17 : void ResetCell::calculate() {
127 :
128 : Pbc & pbc(modifyGlobalPbc());
129 :
130 17 : Tensor box=pbc.getBox();
131 :
132 : // moduli of lattice vectors
133 17 : double a=modulo(box.getRow(0));
134 17 : double b=modulo(box.getRow(1));
135 17 : double c=modulo(box.getRow(2));
136 : // cos-angle between lattice vectors
137 17 : double ab=dotProduct(box.getRow(0),box.getRow(1))/(a*b);
138 17 : double ac=dotProduct(box.getRow(0),box.getRow(2))/(a*c);
139 17 : double bc=dotProduct(box.getRow(1),box.getRow(2))/(b*c);
140 :
141 : // generate a new set of lattice vectors as a lower triangular matrix
142 17 : newbox[0][0]=a;
143 17 : newbox[1][0]=b*ab;
144 17 : newbox[1][1]=std::sqrt(b*b-newbox[1][0]*newbox[1][0]);
145 17 : newbox[2][0]=c*ac;
146 17 : newbox[2][1]=c*(bc-ac*ab)/std::sqrt(1-ab*ab);
147 17 : newbox[2][2]=std::sqrt(c*c-newbox[2][0]*newbox[2][0]-newbox[2][1]*newbox[2][1]);
148 :
149 17 : if(determinant(newbox)*determinant(box)<0) newbox[2][2]=-newbox[2][2];
150 :
151 : // rotation matrix from old to new coordinates
152 17 : rotation=transpose(matmul(inverse(box),newbox));
153 :
154 : // rotate all coordinates
155 3205 : for(unsigned i=0; i<getTotAtoms(); i++) {
156 : Vector & ato (modifyPosition(AtomNumber::index(i)));
157 1594 : ato=matmul(rotation,ato);
158 : }
159 : // rotate box
160 17 : pbc.setBox(newbox);
161 17 : }
162 :
163 17 : void ResetCell::apply() {
164 : // rotate back forces
165 3205 : for(unsigned i=0; i<getTotAtoms(); i++) {
166 : Vector & f(modifyGlobalForce(AtomNumber::index(i)));
167 1594 : f=matmul(transpose(rotation),f);
168 : }
169 :
170 : Tensor& virial(modifyGlobalVirial());
171 : // I have no mathematical derivation for this.
172 : // The reasoning is the following.
173 : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
174 : // The final virial should be rotationally invariant, that is symmetric.
175 : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
176 : // be changed so as to enforce rotational invariance. Thus we here have to
177 : // make the virial matrix symmetric.
178 : // Since h^T is upper triangular, it can be shown that any change in these elements
179 : // will only affect the corresponding elements of the virial matrix.
180 : // Thus, the only possibility is to set the corresponding elements
181 : // of the virial matrix equal to their symmetric ones.
182 : // GB
183 17 : virial[0][1]=virial[1][0];
184 17 : virial[0][2]=virial[2][0];
185 17 : virial[1][2]=virial[2][1];
186 : // rotate back virial
187 17 : virial=matmul(transpose(rotation),matmul(virial,rotation));
188 :
189 :
190 :
191 17 : }
192 :
193 : }
194 4839 : }
|
