Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "tools/PDB.h"
26 : #include "core/PlumedMain.h"
27 :
28 : namespace PLMD {
29 : namespace colvar {
30 :
31 : class RMSDShortcut : public ActionShortcut {
32 : public:
33 : static void registerKeywords(Keywords& keys);
34 : explicit RMSDShortcut(const ActionOptions&);
35 : };
36 :
37 : PLUMED_REGISTER_ACTION(RMSDShortcut,"RMSD")
38 :
39 185 : void RMSDShortcut::registerKeywords(Keywords& keys) {
40 185 : ActionShortcut::registerKeywords( keys );
41 185 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV");
42 185 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
43 185 : keys.add("compulsory","ALIGN","1.0","the weights to use when aligning to the reference structure");
44 185 : keys.add("compulsory","DISPLACE","1.0","the weights to use when calculating the displacement from the reference structure");
45 185 : keys.addFlag("SQUARED",false," This should be setted if you want MSD instead of RMSD ");
46 185 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
47 185 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
48 185 : keys.addFlag("DISPLACEMENT",false,"Calculate the vector of displacements instead of the length of this vector");
49 370 : keys.addInputKeyword("optional","ARG","vector/matrix","instead of using the REFERENCE option you can use this action to specify the labels of two actions that you are calculating the RMSD between");
50 185 : keys.add("compulsory","NUMBER","0","if there are multiple structures in the pdb file you can specify that you want the RMSD from a specific structure by specifying its place in the file here. If NUMBER=0 then the RMSD from all structures are computed");
51 370 : keys.addOutputComponent("disp","DISPLACEMENT","vector/matrix","the vector of displacements for the atoms");
52 370 : keys.addOutputComponent("dist","DISPLACEMENT","scalar/vector","the RMSD distance the atoms have moved");
53 370 : keys.setValueDescription("scalar/vector","the RMSD distance between the instaneous structure and the reference structure/s that were input");
54 185 : keys.addActionNameSuffix("_SCALAR");
55 185 : keys.addActionNameSuffix("_VECTOR");
56 185 : keys.needsAction("PDB2CONSTANT");
57 185 : keys.needsAction("WHOLEMOLECULES");
58 185 : keys.needsAction("POSITION");
59 185 : keys.needsAction("CONCATENATE");
60 185 : keys.needsAction("RMSD");
61 185 : keys.addDOI("10.1107/S0108767388010128");
62 185 : }
63 :
64 147 : RMSDShortcut::RMSDShortcut(const ActionOptions& ao):
65 : Action(ao),
66 147 : ActionShortcut(ao) {
67 : std::string argn;
68 294 : parse("ARG",argn);
69 147 : if( argn.length()>0 ) {
70 82 : readInputLine( getShortcutLabel() + ": RMSD_VECTOR ARG=" + argn + " " + convertInputLineToString() );
71 : return;
72 : }
73 : bool disp;
74 214 : parseFlag("DISPLACEMENT",disp);
75 : std::string reference;
76 108 : parse("REFERENCE",reference);
77 : // Read the reference pdb file
78 106 : PDB pdb;
79 106 : if( !pdb.read(reference,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength()) ) {
80 3 : plumed_merror("missing file " + reference );
81 : }
82 : unsigned frame;
83 105 : parse("NUMBER",frame);
84 : unsigned nf=1;
85 105 : if( frame==0 ) {
86 96 : FILE* fp=std::fopen(reference.c_str(),"r");
87 : bool do_read=true;
88 : nf=0;
89 606 : while ( do_read ) {
90 539 : PDB mypdb;
91 539 : do_read=mypdb.readFromFilepointer(fp,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
92 539 : if( !do_read && nf>0 ) {
93 : break ;
94 : }
95 510 : nf++;
96 539 : }
97 : }
98 : bool nopbc;
99 105 : parseFlag("NOPBC",nopbc);
100 : // Now create the RMSD object
101 105 : std::string rmsd_line = getShortcutLabel() + ": ";
102 105 : if( nf==1 && !disp ) {
103 82 : rmsd_line += "RMSD_SCALAR REFERENCE=" + reference;
104 82 : if(nopbc) {
105 : rmsd_line += " NOPBC";
106 : }
107 : } else {
108 : std::string ffnum;
109 23 : Tools::convert( frame, ffnum );
110 46 : readInputLine( getShortcutLabel() + "_ref: PDB2CONSTANT REFERENCE=" + reference + " NUMBER=" + ffnum );
111 23 : std::vector<AtomNumber> anum( pdb.getAtomNumbers() );
112 23 : if( !nopbc ) {
113 : std::string num;
114 21 : Tools::convert( anum[0].serial(), num );
115 21 : std::string wm_line = "WHOLEMOLECULES ENTITY0=" + num;
116 209 : for(unsigned i=1; i<anum.size(); ++i) {
117 188 : Tools::convert( anum[i].serial(), num );
118 376 : wm_line += "," + num;
119 : }
120 21 : readInputLine( wm_line );
121 : }
122 : std::string num;
123 23 : Tools::convert( anum[0].serial(), num );
124 23 : std::string pos_line = getShortcutLabel() + "_cpos: POSITION NOPBC ATOMS=" + num;
125 223 : for(unsigned i=1; i<anum.size(); ++i) {
126 200 : Tools::convert( anum[i].serial(), num );
127 400 : pos_line += "," + num;
128 : }
129 23 : readInputLine( pos_line );
130 : // Concatenate the three positions together
131 46 : readInputLine( getShortcutLabel() + "_pos: CONCATENATE ARG=" + getShortcutLabel() + "_cpos.x," + getShortcutLabel() + "_cpos.y," + getShortcutLabel() + "_cpos.z");
132 46 : rmsd_line += "RMSD ARG=" + getShortcutLabel() + "_pos," + getShortcutLabel() + "_ref";
133 23 : if( disp ) {
134 : rmsd_line += " DISPLACEMENT";
135 : }
136 : // Now align
137 23 : std::vector<double> align( pdb.getOccupancy() );
138 23 : Tools::convert( align[0], num );
139 23 : rmsd_line += " ALIGN=" + num;
140 223 : for(unsigned i=1; i<align.size(); ++i) {
141 200 : Tools::convert( align[i], num );
142 400 : rmsd_line += "," + num;
143 : }
144 : // And displace
145 23 : std::vector<double> displace( pdb.getBeta() );
146 23 : Tools::convert( displace[0], num );
147 23 : rmsd_line += " DISPLACE=" + num;
148 223 : for(unsigned i=1; i<displace.size(); ++i) {
149 200 : Tools::convert( displace[i], num );
150 400 : rmsd_line += "," + num;
151 : }
152 : }
153 : // And create the RMSD object
154 : bool numder;
155 105 : parseFlag("NUMERICAL_DERIVATIVES",numder);
156 105 : if(numder && nf==1 && !disp ) {
157 : rmsd_line += " NUMERICAL_DERIVATIVES";
158 100 : } else if( numder ) {
159 0 : error("can only use NUMERICAL_DERIVATIVES flag when RMSD is calculating a single scalar value");
160 : }
161 : bool squared;
162 106 : parseFlag("SQUARED",squared);
163 105 : if(squared) {
164 : rmsd_line += " SQUARED";
165 : }
166 : std::string tt;
167 105 : parse("TYPE",tt);
168 211 : readInputLine( rmsd_line + " TYPE=" + tt );
169 110 : }
170 :
171 : }
172 : }
|
