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 "adjmat/AdjacencyMatrixBase.h"
23 : #include "multicolvar/AtomValuePack.h"
24 : #include "HBPammObject.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/KernelFunctions.h"
27 : #include "tools/IFile.h"
28 :
29 : //+PLUMEDOC MATRIX HBPAMM_MATRIX
30 : /*
31 : Adjacency matrix in which two electronegative atoms are adjacent if they are hydrogen bonded
32 :
33 : \par Examples
34 :
35 : */
36 : //+ENDPLUMEDOC
37 :
38 :
39 : namespace PLMD {
40 : namespace pamm {
41 :
42 : class HBPammMatrix : public adjmat::AdjacencyMatrixBase {
43 : private:
44 : unsigned ndonor_types;
45 : double regulariser;
46 : Matrix<HBPammObject> myhb_objs;
47 : public:
48 : /// Create manual
49 : static void registerKeywords( Keywords& keys );
50 : /// Constructor
51 : explicit HBPammMatrix(const ActionOptions&);
52 : /// Setup the connector -- i.e. read in the clusters file
53 : void setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc );
54 : ///
55 : double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
56 : ///
57 : /// Used to check for connections between atoms
58 : bool checkForConnection( const std::vector<double>& myvals ) const {
59 : return !(myvals[1]>epsilon);
60 : }
61 : };
62 :
63 13789 : PLUMED_REGISTER_ACTION(HBPammMatrix,"HBPAMM_MATRIX")
64 :
65 6 : void HBPammMatrix::registerKeywords( Keywords& keys ) {
66 6 : adjmat::AdjacencyMatrixBase::registerKeywords( keys );
67 12 : keys.add("atoms-1","SITES","The list of atoms which can be part of a hydrogen bond. When this command is used the set of atoms that can donate a "
68 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified "
69 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
70 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
71 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
72 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
73 12 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
74 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
75 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
76 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
77 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
78 12 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
79 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
80 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
81 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
82 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
83 12 : keys.add("atoms","HYDROGENS","The list of hydrogen atoms that can form part of a hydrogen bond. The atoms must be specified using a comma separated list, "
84 : "an index range or by using a \\ref GROUP");
85 12 : keys.add("numbered","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
86 12 : keys.reset_style("CLUSTERS","compulsory");
87 6 : keys.use("SUM");
88 12 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
89 6 : }
90 :
91 :
92 2 : HBPammMatrix::HBPammMatrix(const ActionOptions& ao):
93 : Action(ao),
94 2 : AdjacencyMatrixBase(ao) {
95 4 : readMaxThreeSpeciesMatrix("SITES", "DONORS", "ACCEPTORS", "HYDROGENS", false );
96 : // Retrieve dimensions of hbonding matrix and resize
97 : unsigned nrows, ncols;
98 2 : retrieveTypeDimensions( nrows, ncols, ndonor_types );
99 2 : myhb_objs.resize( nrows, ncols );
100 : // Read in the regularisation parameter
101 2 : parse("REGULARISE",regulariser);
102 : // Read in the switching functions
103 2 : parseConnectionDescriptions("CLUSTERS",false,ndonor_types);
104 :
105 : // Find cutoff for link cells
106 2 : double sfmax=0;
107 4 : for(unsigned i=0; i<myhb_objs.ncols(); ++i) {
108 4 : for(unsigned j=i; j<myhb_objs.nrows(); ++j) {
109 2 : double rcut=myhb_objs(i,j).get_cutoff();
110 2 : if( rcut>sfmax ) {
111 2 : sfmax=rcut;
112 : }
113 : }
114 : }
115 2 : setLinkCellCutoff( sfmax );
116 2 : }
117 :
118 2 : void HBPammMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
119 2 : log.printf(" reading definition of hydrogen bond between between type %u and %u from file %s \n",i,j,desc[0].c_str() );
120 2 : plumed_assert( desc.size()==1 );
121 : std::string errors;
122 2 : if( i==j ) {
123 2 : myhb_objs( i, j ).setup( desc[0], regulariser, this, errors );
124 : } else {
125 0 : myhb_objs( i, j ).setup( desc[0], regulariser, this, errors );
126 0 : myhb_objs( j, i ).setup( desc[0], regulariser, this, errors );
127 : }
128 2 : if( errors.length()>0 ) {
129 0 : error( errors );
130 : }
131 2 : }
132 :
133 4038 : double HBPammMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
134 4038 : Vector d_da = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
135 4038 : double md_da = d_da.modulo(); // acceptor - donor
136 :
137 : // Get the base colvar numbers
138 : unsigned ano, dno = getBaseColvarNumber( myatoms.getIndex(0) );
139 4038 : if( ndonor_types==0 ) {
140 : ano = getBaseColvarNumber( myatoms.getIndex(1) );
141 : } else {
142 0 : ano = getBaseColvarNumber( myatoms.getIndex(1) ) - ndonor_types;
143 : }
144 :
145 : double value=0;
146 4038 : if( myatoms.getNumberOfAtoms()>3 ) {
147 : const HBPammObject& myhb=myhb_objs(dno,ano);
148 520170 : for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) {
149 516132 : value+=myhb.evaluate( 0, 1, i, d_da, md_da, myatoms );
150 : }
151 : } else {
152 : plumed_dbg_assert( myatoms.getNumberOfAtoms()==3 );
153 0 : value=myhb_objs(dno,ano).evaluate( 0, 1, 2, d_da, md_da, myatoms );
154 : }
155 :
156 4038 : return value;
157 : }
158 :
159 : }
160 : }
|
