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 220 : void RMSDShortcut::registerKeywords(Keywords& keys) {
40 220 : ActionShortcut::registerKeywords( keys );
41 220 : keys.add("compulsory","REFERENCE","a file in pdb format containing the reference structure and the atoms involved in the CV");
42 220 : keys.add("compulsory","TYPE","SIMPLE","the manner in which RMSD alignment is performed. Should be OPTIMAL or SIMPLE.");
43 220 : keys.add("compulsory","ALIGN","1.0","the weights to use when aligning to the reference structure");
44 220 : keys.add("compulsory","DISPLACE","1.0","the weights to use when calculating the displacement from the reference structure");
45 220 : keys.addFlag("SQUARED",false," This should be set if you want MSD instead of RMSD ");
46 220 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
47 220 : keys.addFlag("NUMERICAL_DERIVATIVES", false, "calculate the derivatives for these quantities numerically");
48 220 : keys.addFlag("DISPLACEMENT",false,"Calculate the vector of displacements instead of the length of this vector");
49 440 : 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 220 : 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 440 : keys.addOutputComponent("disp","DISPLACEMENT","vector/matrix","the vector of displacements for the atoms");
52 440 : keys.addOutputComponent("dist","DISPLACEMENT","scalar/vector","the RMSD distance the atoms have moved");
53 440 : keys.setValueDescription("scalar/vector","the RMSD distance between the instantaneous structure and the reference structure/s that were input");
54 220 : keys.addActionNameSuffix("_SCALAR");
55 220 : keys.addActionNameSuffix("_VECTOR");
56 220 : keys.needsAction("PDB2CONSTANT");
57 220 : keys.needsAction("WHOLEMOLECULES");
58 220 : keys.needsAction("POSITION");
59 220 : keys.needsAction("CONCATENATE");
60 220 : keys.needsAction("RMSD");
61 220 : keys.addDOI("10.1107/S0108767388010128");
62 220 : }
63 :
64 169 : RMSDShortcut::RMSDShortcut(const ActionOptions& ao):
65 : Action(ao),
66 169 : ActionShortcut(ao) {
67 : std::string argn;
68 338 : parse("ARG",argn);
69 169 : if( argn.length()>0 ) {
70 165 : readInputLine( getShortcutLabel() + ": RMSD_VECTOR ARG=" + argn + " " + convertInputLineToString() );
71 : return;
72 : }
73 : bool disp;
74 231 : parseFlag("DISPLACEMENT",disp);
75 : std::string reference;
76 117 : parse("REFERENCE",reference);
77 : // Read the reference pdb file
78 114 : PDB pdb;
79 114 : if( !pdb.read(reference,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength()) ) {
80 3 : plumed_merror("missing file " + reference );
81 : }
82 : unsigned frame;
83 112 : parse("NUMBER",frame);
84 : unsigned nf=1;
85 112 : if( frame==0 ) {
86 103 : FILE* fp=std::fopen(reference.c_str(),"r");
87 103 : if(!fp) {
88 0 : plumed_merror("cannot open file " + reference );
89 : }
90 : bool do_read=true;
91 : nf=0;
92 761 : while ( do_read ) {
93 694 : PDB mypdb;
94 694 : do_read=mypdb.readFromFilepointer(fp,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
95 694 : if( !do_read && nf>0 ) {
96 : break ;
97 : }
98 658 : nf++;
99 694 : }
100 103 : std::fclose(fp);
101 : }
102 : bool nopbc;
103 112 : parseFlag("NOPBC",nopbc);
104 : // Now create the RMSD object
105 112 : std::string rmsd_line = getShortcutLabel() + ": ";
106 112 : if( nf==1 && !disp ) {
107 84 : rmsd_line += "RMSD_SCALAR REFERENCE=" + reference;
108 83 : if(nopbc) {
109 : rmsd_line += " NOPBC";
110 : }
111 : } else {
112 : std::string ffnum;
113 29 : Tools::convert( frame, ffnum );
114 58 : readInputLine( getShortcutLabel() + "_ref: PDB2CONSTANT REFERENCE=" + reference + " NUMBER=" + ffnum );
115 29 : std::vector<AtomNumber> anum( pdb.getAtomNumbers() );
116 29 : if( !nopbc ) {
117 : std::string num;
118 27 : Tools::convert( anum[0].serial(), num );
119 27 : std::string wm_line = "WHOLEMOLECULES ENTITY0=" + num;
120 286 : for(unsigned i=1; i<anum.size(); ++i) {
121 259 : Tools::convert( anum[i].serial(), num );
122 518 : wm_line += "," + num;
123 : }
124 27 : readInputLine( wm_line );
125 : }
126 : std::string num;
127 29 : Tools::convert( anum[0].serial(), num );
128 29 : std::string pos_line = getShortcutLabel() + "_cpos: POSITION NOPBC ATOMS=" + num;
129 300 : for(unsigned i=1; i<anum.size(); ++i) {
130 271 : Tools::convert( anum[i].serial(), num );
131 542 : pos_line += "," + num;
132 : }
133 29 : readInputLine( pos_line );
134 : // Concatenate the three positions together
135 116 : readInputLine( getShortcutLabel() + "_pos: CONCATENATE ARG=" + getShortcutLabel() + "_cpos.x," + getShortcutLabel() + "_cpos.y," + getShortcutLabel() + "_cpos.z");
136 58 : rmsd_line += "RMSD ARG=" + getShortcutLabel() + "_pos," + getShortcutLabel() + "_ref";
137 29 : if( disp ) {
138 : rmsd_line += " DISPLACEMENT";
139 : }
140 : // Now align
141 29 : std::vector<double> align( pdb.getOccupancy() );
142 29 : Tools::convert( align[0], num );
143 29 : rmsd_line += " ALIGN=" + num;
144 300 : for(unsigned i=1; i<align.size(); ++i) {
145 271 : Tools::convert( align[i], num );
146 542 : rmsd_line += "," + num;
147 : }
148 : // And displace
149 29 : std::vector<double> displace( pdb.getBeta() );
150 29 : Tools::convert( displace[0], num );
151 29 : rmsd_line += " DISPLACE=" + num;
152 300 : for(unsigned i=1; i<displace.size(); ++i) {
153 271 : Tools::convert( displace[i], num );
154 542 : rmsd_line += "," + num;
155 : }
156 58 : }
157 : // And create the RMSD object
158 : bool numder;
159 112 : parseFlag("NUMERICAL_DERIVATIVES",numder);
160 112 : if(numder && nf==1 && !disp ) {
161 : rmsd_line += " NUMERICAL_DERIVATIVES";
162 107 : } else if( numder ) {
163 0 : error("can only use NUMERICAL_DERIVATIVES flag when RMSD is calculating a single scalar value");
164 : }
165 : bool squared;
166 113 : parseFlag("SQUARED",squared);
167 112 : if(squared) {
168 : rmsd_line += " SQUARED";
169 : }
170 : std::string tt;
171 112 : parse("TYPE",tt);
172 224 : readInputLine( rmsd_line + " TYPE=" + tt );
173 203 : }
174 :
175 : }
176 : }
|
