Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 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 "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "tools/Vector.h"
28 : #include "tools/Matrix.h"
29 : #include "tools/AtomNumber.h"
30 : #include "tools/Tools.h"
31 : #include "core/PbcAction.h"
32 : #include "tools/Pbc.h"
33 :
34 : namespace PLMD {
35 : namespace generic {
36 :
37 : //+PLUMEDOC GENERIC RESET_CELL
38 : /*
39 : This action is used to rotate the full cell
40 :
41 : This can be used to modify the periodic box. Notice that
42 : this is done at fixed scaled coordinates,
43 : so that also atomic coordinates for the entire system are affected.
44 : To see what effect try
45 : the \ref DUMPATOMS directive to output the atomic positions.
46 :
47 : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
48 : after rotation. See also \ref FIT_TO_TEMPLATE
49 :
50 : Currently, only TYPE=TRIANGULAR is implemented, which allows one to reset
51 : the cell to a lower triangular one. Namely, a proper rotation is found that allows
52 : rotating the box so that the first lattice vector is in the form (ax,0,0),
53 : the second lattice vector is in the form (bx,by,0), and the third lattice vector is
54 : arbitrary.
55 :
56 : \attention
57 : The implementation of this action is available but should be considered in testing phase. Please report any
58 : strange behavior.
59 :
60 : \attention
61 : This directive modifies the stored position at the precise moment
62 : it is executed. This means that only collective variables
63 : which are below it in the input script will see the corrected positions.
64 : Unless you
65 : know exactly what you are doing, leave the default stride (1), so that
66 : this action is performed at every MD step.
67 :
68 : \par Examples
69 :
70 : Reset cell to be triangular after a rototranslational fit
71 : \plumedfile
72 : DUMPATOMS FILE=dump-original.xyz ATOMS=1-20
73 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=OPTIMAL
74 : DUMPATOMS FILE=dump-fit.xyz ATOMS=1-20
75 : RESET_CELL TYPE=TRIANGULAR
76 : DUMPATOMS FILE=dump-reset.xyz ATOMS=1-20
77 : \endplumedfile
78 :
79 : The reference file for the FIT_TO_TEMPLATE is just a normal pdb file with the format shown below:
80 :
81 : \auxfile{ref.pdb}
82 : ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H
83 : ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C
84 : ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H
85 : ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H
86 : ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O
87 : END
88 : \endauxfile
89 :
90 : */
91 : //+ENDPLUMEDOC
92 :
93 :
94 : class ResetCell:
95 : public ActionPilot,
96 : public ActionAtomistic {
97 : std::string type;
98 : Tensor rotation,newbox;
99 : Value* boxValue;
100 : PbcAction* pbc_action;
101 : public:
102 : explicit ResetCell(const ActionOptions&ao);
103 : static void registerKeywords( Keywords& keys );
104 : void calculate() override;
105 : void apply() override;
106 : };
107 :
108 : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
109 :
110 6 : void ResetCell::registerKeywords( Keywords& keys ) {
111 6 : Action::registerKeywords( keys );
112 6 : ActionAtomistic::registerKeywords( keys );
113 12 : 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!");
114 12 : keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
115 6 : }
116 :
117 2 : ResetCell::ResetCell(const ActionOptions&ao):
118 : Action(ao),
119 : ActionPilot(ao),
120 2 : ActionAtomistic(ao) {
121 2 : type.assign("TRIANGULAR");
122 2 : parse("TYPE",type);
123 :
124 2 : log<<" type: "<<type<<"\n";
125 2 : if(type!="TRIANGULAR") {
126 0 : error("undefined type "+type);
127 : }
128 :
129 2 : pbc_action=plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
130 2 : if( !pbc_action ) {
131 0 : error("cannot reset cell if box has not been set");
132 : }
133 2 : boxValue=pbc_action->copyOutput(0);
134 2 : }
135 :
136 :
137 17 : void ResetCell::calculate() {
138 17 : Pbc & pbc(pbc_action->getPbc());
139 17 : Tensor box=pbc.getBox();
140 :
141 : // moduli of lattice vectors
142 17 : double a=modulo(box.getRow(0));
143 17 : double b=modulo(box.getRow(1));
144 17 : double c=modulo(box.getRow(2));
145 : // cos-angle between lattice vectors
146 17 : double ab=dotProduct(box.getRow(0),box.getRow(1))/(a*b);
147 17 : double ac=dotProduct(box.getRow(0),box.getRow(2))/(a*c);
148 17 : double bc=dotProduct(box.getRow(1),box.getRow(2))/(b*c);
149 :
150 : // generate a new set of lattice vectors as a lower triangular matrix
151 17 : newbox[0][0]=a;
152 17 : newbox[1][0]=b*ab;
153 17 : newbox[1][1]=std::sqrt(b*b-newbox[1][0]*newbox[1][0]);
154 17 : newbox[2][0]=c*ac;
155 17 : newbox[2][1]=c*(bc-ac*ab)/std::sqrt(1-ab*ab);
156 17 : newbox[2][2]=std::sqrt(c*c-newbox[2][0]*newbox[2][0]-newbox[2][1]*newbox[2][1]);
157 :
158 17 : if(determinant(newbox)*determinant(box)<0) {
159 1 : newbox[2][2]=-newbox[2][2];
160 : }
161 :
162 : // rotation matrix from old to new coordinates
163 17 : rotation=transpose(matmul(inverse(box),newbox));
164 :
165 : // rotate all coordinates
166 17 : unsigned nat = getTotAtoms();
167 1623 : for(unsigned i=0; i<nat; i++) {
168 1606 : std::pair<std::size_t,std::size_t> a = getValueIndices( AtomNumber::index(i));
169 1606 : Vector ato=matmul(rotation,getGlobalPosition(a));
170 1606 : setGlobalPosition(a,ato);
171 : }
172 : // rotate box
173 17 : pbc.setBox(newbox);
174 17 : }
175 :
176 17 : void ResetCell::apply() {
177 : // rotate back forces
178 17 : unsigned nat = getTotAtoms();
179 1623 : for(unsigned i=0; i<nat; i++) {
180 1606 : std::pair<std::size_t,std::size_t> a = getValueIndices( AtomNumber::index(i));
181 1606 : Vector f=getForce(a);
182 1606 : Vector nf=matmul(transpose(rotation),f);
183 1606 : addForce(a, nf-f );
184 : }
185 :
186 : // I have no mathematical derivation for this.
187 : // The reasoning is the following.
188 : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
189 : // The final virial should be rotationally invariant, that is symmetric.
190 : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
191 : // be changed so as to enforce rotational invariance. Thus we here have to
192 : // make the virial matrix symmetric.
193 : // Since h^T is upper triangular, it can be shown that any change in these elements
194 : // will only affect the corresponding elements of the virial matrix.
195 : // Thus, the only possibility is to set the corresponding elements
196 : // of the virial matrix equal to their symmetric ones.
197 : // GB
198 17 : Tensor virial;
199 68 : for(unsigned i=0; i<3; ++i)
200 204 : for(unsigned j=0; j<3; ++j) {
201 153 : virial[i][j]=boxValue->getForce( 3*i+j );
202 : }
203 17 : virial[0][1]=virial[1][0];
204 17 : virial[0][2]=virial[2][0];
205 17 : virial[1][2]=virial[2][1];
206 : // rotate back virial
207 17 : virial=matmul(transpose(rotation),matmul(virial,rotation));
208 17 : boxValue->clearInputForce();
209 68 : for(unsigned i=0; i<3; ++i)
210 204 : for(unsigned j=0; j<3; ++j) {
211 153 : boxValue->addForce( 3*i+j, virial(i,j) );
212 : }
213 :
214 :
215 17 : }
216 :
217 : }
218 : }
|
