Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "AlignedMatrixBase.h"
23 : #include "multicolvar/AtomValuePack.h"
24 : #include "tools/SwitchingFunction.h"
25 : #include "tools/Matrix.h"
26 :
27 : namespace PLMD {
28 : namespace adjmat {
29 :
30 9 : void AlignedMatrixBase::registerKeywords( Keywords& keys ) {
31 9 : AdjacencyMatrixBase::registerKeywords( keys );
32 18 : keys.add("atoms","ATOMS","The list of molecules for which you would like to calculate the contact matrix. The molecules involved must "
33 : "have an orientation so your list will be a list of the labels of \\ref mcolv or \\ref multicolvarfunction "
34 : "as PLUMED calculates the orientations of molecules within these operations. Please note also that the majority "
35 : "of \\ref mcolv and \\ref multicolvarfunction do not calculate a molecular orientation.");
36 18 : keys.add("atoms-1","ATOMSA","The list of molecules that you would like to use for the rows of the contact matrix. The molecules involved must "
37 : "have an orientation so your list will be a list of the labels of \\ref mcolv or \\ref multicolvarfunction "
38 : "as PLUMED calculates the orientations of molecules within these operations. Please note also that the majority "
39 : "of \\ref mcolv and \\ref multicolvarfunction do not calculate a molecular orientation.");
40 18 : keys.add("atoms-1","ATOMSB","The list of molecules that you would like to use for the columns of the contact matrix. The molecules involved must "
41 : "have an orientation so your list will be a list of the labels of \\ref mcolv or \\ref multicolvarfunction "
42 : "as PLUMED calculates the orientations of molecules within these operations. Please note also that the majority "
43 : "of \\ref mcolv and \\ref multicolvarfunction do not calculate a molecular orientation.");
44 18 : keys.add("numbered","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
45 : "The following provides information on the \\ref switchingfunction that are available.");
46 9 : }
47 :
48 1 : AlignedMatrixBase::AlignedMatrixBase( const ActionOptions& ao ):
49 : Action(ao),
50 1 : AdjacencyMatrixBase(ao) {
51 : // Read in the atomic positions
52 2 : readMaxTwoSpeciesMatrix( "ATOMS", "ATOMSA", "ATOMSB", true );
53 : unsigned nrows, ncols;
54 1 : retrieveTypeDimensions( nrows, ncols, ncol_t );
55 1 : if( mybasemulticolvars.size()==0 ) {
56 0 : error("cannot use atom indices as input to this variable / input was not specified");
57 : }
58 1 : if( getSizeOfInputVectors()<3 ) {
59 0 : error("base multicolvars do not calculate an orientation");
60 : }
61 : // Read in the switching function
62 1 : switchingFunction.resize( nrows, ncols );
63 2 : parseConnectionDescriptions("SWITCH",false,ncol_t);
64 :
65 : // Find the largest sf cutoff
66 1 : double sfmax=switchingFunction(0,0).get_dmax();
67 2 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
68 2 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
69 1 : double tsf=switchingFunction(i,j).get_dmax();
70 1 : if( tsf>sfmax ) {
71 0 : sfmax=tsf;
72 : }
73 : }
74 : }
75 : // And set the link cell cutoff
76 1 : setLinkCellCutoff( sfmax );
77 1 : }
78 :
79 2 : void AlignedMatrixBase::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
80 2 : plumed_assert( id<2 );
81 2 : if( id==0 ) {
82 1 : plumed_assert( desc.size()==1 );
83 : std::string errors;
84 1 : switchingFunction(j,i).set(desc[0],errors);
85 1 : if( errors.length()!=0 ) {
86 0 : error("problem reading switching function in SWITCH keywrd description " + errors);
87 : }
88 1 : if( j!=i) {
89 0 : switchingFunction(i,j).set(desc[0],errors);
90 : }
91 2 : log.printf(" %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
92 1 : } else if( id==1 ) {
93 1 : readOrientationConnector( i, j, desc );
94 : }
95 2 : }
96 :
97 57063 : double AlignedMatrixBase::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
98 57063 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
99 171189 : if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).get_dmax2() ) {
100 5806 : return 1.0;
101 : }
102 : return 0.0;
103 : }
104 :
105 5806 : double AlignedMatrixBase::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
106 5806 : unsigned ncomp=getSizeOfInputVectors();
107 5806 : Vector ddistance;
108 5806 : std::vector<double> orient0(ncomp), orient1(ncomp), dorient0(ncomp), dorient1(ncomp);
109 5806 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
110 5806 : getInputData( 0, true, myatoms, orient0 );
111 5806 : getInputData( 1, true, myatoms, orient1 );
112 17418 : double f_dot = computeVectorFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t,
113 : distance, orient0, orient1, ddistance, dorient0, dorient1 );
114 :
115 : // Retrieve the weight of the connection
116 11612 : double dfunc, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).calculate( distance.modulo(), dfunc );
117 :
118 5806 : if( !doNotCalculateDerivatives() ) {
119 0 : addAtomDerivatives( 1, 0, (-dfunc)*f_dot*distance - sw*ddistance, myatoms );
120 0 : addAtomDerivatives( 1, 1, (+dfunc)*f_dot*distance + sw*ddistance, myatoms );
121 0 : myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct( distance, ddistance ) );
122 :
123 : // Add derivatives of orientation
124 0 : for(unsigned k=2; k<orient0.size(); ++k) {
125 0 : dorient0[k]*=sw;
126 0 : dorient1[k]*=sw;
127 : }
128 0 : mergeInputDerivatives( 1, 2, orient0.size(), 0, dorient0, getInputDerivatives( 0, true, myatoms ), myatoms );
129 0 : mergeInputDerivatives( 1, 2, orient1.size(), 1, dorient1, getInputDerivatives( 1, true, myatoms ), myatoms );
130 : }
131 11612 : return sw*f_dot;
132 : }
133 :
134 : }
135 : }
136 :
|
