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 { return !(myvals[1]>epsilon); }
59 : };
60 :
61 10423 : PLUMED_REGISTER_ACTION(HBPammMatrix,"HBPAMM_MATRIX")
62 :
63 3 : void HBPammMatrix::registerKeywords( Keywords& keys ) {
64 3 : adjmat::AdjacencyMatrixBase::registerKeywords( keys );
65 6 : 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 "
66 : "hydrogen bond is assumed to be the same as the set of atoms that can form hydrogen bonds. The atoms involved must be specified "
67 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
68 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
69 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
70 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
71 6 : keys.add("atoms-2","DONORS","The list of atoms which can donate a hydrogen bond. The atoms involved must be specified "
72 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
73 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
74 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
75 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
76 6 : keys.add("atoms-2","ACCEPTORS","The list of atoms which can accept a hydrogen bond. The atoms involved must be specified "
77 : "as a list of labels of \\ref mcolv or labels of a \\ref multicolvarfunction actions. If you would just like to use "
78 : "the atomic positions you can use a \\ref DENSITY command to specify a group of atoms. Specifying your atomic positions using labels of "
79 : "other \\ref mcolv or \\ref multicolvarfunction commands is useful, however, as you can then exploit a much wider "
80 : "variety of functions of the contact matrix as described in \\ref contactmatrix");
81 6 : 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, "
82 : "an index range or by using a \\ref GROUP");
83 6 : keys.add("numbered","CLUSTERS","the name of the file that contains the definitions of all the kernels for PAMM");
84 6 : keys.reset_style("CLUSTERS","compulsory"); keys.use("SUM");
85 6 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
86 3 : }
87 :
88 :
89 2 : HBPammMatrix::HBPammMatrix(const ActionOptions& ao):
90 : Action(ao),
91 2 : AdjacencyMatrixBase(ao)
92 : {
93 4 : readMaxThreeSpeciesMatrix("SITES", "DONORS", "ACCEPTORS", "HYDROGENS", false );
94 : // Retrieve dimensions of hbonding matrix and resize
95 2 : unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ndonor_types );
96 2 : myhb_objs.resize( nrows, ncols );
97 : // Read in the regularisation parameter
98 2 : parse("REGULARISE",regulariser);
99 : // Read in the switching functions
100 2 : parseConnectionDescriptions("CLUSTERS",false,ndonor_types);
101 :
102 : // Find cutoff for link cells
103 2 : double sfmax=0;
104 4 : for(unsigned i=0; i<myhb_objs.ncols(); ++i) {
105 4 : for(unsigned j=i; j<myhb_objs.nrows(); ++j) {
106 2 : double rcut=myhb_objs(i,j).get_cutoff();
107 2 : if( rcut>sfmax ) { sfmax=rcut; }
108 : }
109 : }
110 2 : setLinkCellCutoff( sfmax );
111 2 : }
112 :
113 2 : void HBPammMatrix::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
114 2 : log.printf(" reading definition of hydrogen bond between between type %u and %u from file %s \n",i,j,desc[0].c_str() );
115 2 : plumed_assert( desc.size()==1 ); std::string errors;
116 2 : if( i==j ) {
117 2 : myhb_objs( i, j ).setup( desc[0], regulariser, this, errors );
118 : } else {
119 0 : myhb_objs( i, j ).setup( desc[0], regulariser, this, errors );
120 0 : myhb_objs( j, i ).setup( desc[0], regulariser, this, errors );
121 : }
122 2 : if( errors.length()>0 ) error( errors );
123 2 : }
124 :
125 4038 : double HBPammMatrix::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
126 4038 : Vector d_da = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) ); double md_da = d_da.modulo(); // acceptor - donor
127 :
128 : // Get the base colvar numbers
129 : unsigned ano, dno = getBaseColvarNumber( myatoms.getIndex(0) );
130 4038 : if( ndonor_types==0 ) ano = getBaseColvarNumber( myatoms.getIndex(1) );
131 0 : else ano = getBaseColvarNumber( myatoms.getIndex(1) ) - ndonor_types;
132 :
133 : double value=0;
134 4038 : if( myatoms.getNumberOfAtoms()>3 ) {
135 : const HBPammObject& myhb=myhb_objs(dno,ano);
136 520170 : for(unsigned i=2; i<myatoms.getNumberOfAtoms(); ++i) {
137 516132 : value+=myhb.evaluate( 0, 1, i, d_da, md_da, myatoms );
138 : }
139 : } else {
140 : plumed_dbg_assert( myatoms.getNumberOfAtoms()==3 );
141 0 : value=myhb_objs(dno,ano).evaluate( 0, 1, 2, d_da, md_da, myatoms );
142 : }
143 :
144 4038 : return value;
145 : }
146 :
147 : }
148 : }
|
