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 : Rotating the full cell is useful if you want 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 using [DUMPATOMS](DUMPATOMS.md) directive to output the atomic positions.
45 :
46 : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
47 : after rotation. You can read the documentation for [FIT_TO_TEMPLATE](FIT_TO_TEMPLATE.md) for more
48 : detail.
49 :
50 : Currently, only TYPE=TRIANGULAR is implemented, which allows one to reset
51 : the cell to a lower triangular one. This command finds a proper rotation
52 : rotates 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 : > [!CAUTION]
57 : > The implementation of this action is available but should be considered in testing phase. Please report any
58 : > strange behavior.
59 :
60 : > [!CAUTION]
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 : ## Examples
69 :
70 : Reset cell to be triangular after a rototranslational fit
71 :
72 : ```plumed
73 : #SETTINGS INPUTFILES=regtest/basic/rt63/align.pdb
74 : DUMPATOMS FILE=dump-original.xyz ATOMS=1-20
75 : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=regtest/basic/rt63/align.pdb TYPE=OPTIMAL
76 : DUMPATOMS FILE=dump-fit.xyz ATOMS=1-20
77 : RESET_CELL TYPE=TRIANGULAR
78 : DUMPATOMS FILE=dump-reset.xyz ATOMS=1-20
79 : ```
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 :
85 : class ResetCell:
86 : public ActionPilot,
87 : public ActionAtomistic {
88 : std::string type;
89 : Tensor rotation,newbox;
90 : Value* boxValue;
91 : PbcAction* pbc_action;
92 : public:
93 : explicit ResetCell(const ActionOptions&ao);
94 : static void registerKeywords( Keywords& keys );
95 : void calculate() override;
96 : void apply() override;
97 : };
98 :
99 : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
100 :
101 4 : void ResetCell::registerKeywords( Keywords& keys ) {
102 4 : Action::registerKeywords( keys );
103 4 : ActionAtomistic::registerKeywords( keys );
104 4 : 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!");
105 4 : keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
106 4 : }
107 :
108 2 : ResetCell::ResetCell(const ActionOptions&ao):
109 : Action(ao),
110 : ActionPilot(ao),
111 2 : ActionAtomistic(ao) {
112 2 : type.assign("TRIANGULAR");
113 2 : parse("TYPE",type);
114 :
115 2 : log<<" type: "<<type<<"\n";
116 2 : if(type!="TRIANGULAR") {
117 0 : error("undefined type "+type);
118 : }
119 :
120 2 : pbc_action=plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
121 2 : if( !pbc_action ) {
122 0 : error("cannot reset cell if box has not been set");
123 : }
124 2 : boxValue=pbc_action->copyOutput(0);
125 2 : }
126 :
127 :
128 17 : void ResetCell::calculate() {
129 17 : Pbc & pbc(pbc_action->getPbc());
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) {
150 1 : newbox[2][2]=-newbox[2][2];
151 : }
152 :
153 : // rotation matrix from old to new coordinates
154 17 : rotation=transpose(matmul(inverse(box),newbox));
155 :
156 : // rotate all coordinates
157 17 : unsigned nat = getTotAtoms();
158 1623 : for(unsigned i=0; i<nat; i++) {
159 1606 : std::pair<std::size_t,std::size_t> a = getValueIndices( AtomNumber::index(i));
160 1606 : Vector ato=matmul(rotation,getGlobalPosition(a));
161 1606 : setGlobalPosition(a,ato);
162 : }
163 : // rotate box
164 17 : pbc.setBox(newbox);
165 17 : }
166 :
167 17 : void ResetCell::apply() {
168 : // rotate back forces
169 17 : unsigned nat = getTotAtoms();
170 1623 : for(unsigned i=0; i<nat; i++) {
171 1606 : std::pair<std::size_t,std::size_t> a = getValueIndices( AtomNumber::index(i));
172 1606 : Vector f=getForce(a);
173 1606 : Vector nf=matmul(transpose(rotation),f);
174 1606 : addForce(a, nf-f );
175 : }
176 :
177 : // I have no mathematical derivation for this.
178 : // The reasoning is the following.
179 : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
180 : // The final virial should be rotationally invariant, that is symmetric.
181 : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
182 : // be changed so as to enforce rotational invariance. Thus we here have to
183 : // make the virial matrix symmetric.
184 : // Since h^T is upper triangular, it can be shown that any change in these elements
185 : // will only affect the corresponding elements of the virial matrix.
186 : // Thus, the only possibility is to set the corresponding elements
187 : // of the virial matrix equal to their symmetric ones.
188 : // GB
189 17 : Tensor virial;
190 68 : for(unsigned i=0; i<3; ++i)
191 204 : for(unsigned j=0; j<3; ++j) {
192 153 : virial[i][j]=boxValue->getForce( 3*i+j );
193 : }
194 17 : virial[0][1]=virial[1][0];
195 17 : virial[0][2]=virial[2][0];
196 17 : virial[1][2]=virial[2][1];
197 : // rotate back virial
198 17 : virial=matmul(transpose(rotation),matmul(virial,rotation));
199 17 : boxValue->clearInputForce();
200 68 : for(unsigned i=0; i<3; ++i)
201 204 : for(unsigned j=0; j<3; ++j) {
202 153 : boxValue->addForce( 3*i+j, virial(i,j) );
203 : }
204 :
205 :
206 17 : }
207 :
208 : }
209 : }
|
