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 3 : void AlignedMatrixBase::registerKeywords( Keywords& keys ) {
31 3 : AdjacencyMatrixBase::registerKeywords( keys );
32 6 : 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 6 : 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 6 : 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 6 : 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 3 : }
47 :
48 1 : AlignedMatrixBase::AlignedMatrixBase( const ActionOptions& ao ):
49 : Action(ao),
50 1 : AdjacencyMatrixBase(ao)
51 : {
52 : // Read in the atomic positions
53 2 : readMaxTwoSpeciesMatrix( "ATOMS", "ATOMSA", "ATOMSB", true );
54 1 : unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ncol_t );
55 1 : if( mybasemulticolvars.size()==0 ) error("cannot use atom indices as input to this variable / input was not specified");
56 1 : if( getSizeOfInputVectors()<3 ) error("base multicolvars do not calculate an orientation");
57 : // Read in the switching function
58 1 : switchingFunction.resize( nrows, ncols );
59 2 : parseConnectionDescriptions("SWITCH",false,ncol_t);
60 :
61 : // Find the largest sf cutoff
62 1 : double sfmax=switchingFunction(0,0).get_dmax();
63 2 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
64 2 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
65 1 : double tsf=switchingFunction(i,j).get_dmax();
66 1 : if( tsf>sfmax ) sfmax=tsf;
67 : }
68 : }
69 : // And set the link cell cutoff
70 1 : setLinkCellCutoff( sfmax );
71 1 : }
72 :
73 2 : void AlignedMatrixBase::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
74 2 : plumed_assert( id<2 );
75 2 : if( id==0 ) {
76 1 : plumed_assert( desc.size()==1 ); std::string errors; switchingFunction(j,i).set(desc[0],errors);
77 1 : if( errors.length()!=0 ) error("problem reading switching function in SWITCH keywrd description " + errors);
78 1 : if( j!=i) switchingFunction(i,j).set(desc[0],errors);
79 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() );
80 1 : } else if( id==1 ) {
81 1 : readOrientationConnector( i, j, desc );
82 : }
83 2 : }
84 :
85 57063 : double AlignedMatrixBase::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
86 57063 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
87 171189 : if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).get_dmax2() ) return 1.0;
88 : return 0.0;
89 : }
90 :
91 5806 : double AlignedMatrixBase::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
92 5806 : unsigned ncomp=getSizeOfInputVectors(); Vector ddistance;
93 5806 : std::vector<double> orient0(ncomp), orient1(ncomp), dorient0(ncomp), dorient1(ncomp);
94 5806 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
95 5806 : getInputData( 0, true, myatoms, orient0 ); getInputData( 1, true, myatoms, orient1 );
96 17418 : double f_dot = computeVectorFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t,
97 : distance, orient0, orient1, ddistance, dorient0, dorient1 );
98 :
99 : // Retrieve the weight of the connection
100 11612 : double dfunc, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).calculate( distance.modulo(), dfunc );
101 :
102 5806 : if( !doNotCalculateDerivatives() ) {
103 0 : addAtomDerivatives( 1, 0, (-dfunc)*f_dot*distance - sw*ddistance, myatoms );
104 0 : addAtomDerivatives( 1, 1, (+dfunc)*f_dot*distance + sw*ddistance, myatoms );
105 0 : myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct( distance, ddistance ) );
106 :
107 : // Add derivatives of orientation
108 0 : for(unsigned k=2; k<orient0.size(); ++k) { dorient0[k]*=sw; dorient1[k]*=sw; }
109 0 : mergeInputDerivatives( 1, 2, orient0.size(), 0, dorient0, getInputDerivatives( 0, true, myatoms ), myatoms );
110 0 : mergeInputDerivatives( 1, 2, orient1.size(), 1, dorient1, getInputDerivatives( 1, true, myatoms ), myatoms );
111 : }
112 11612 : return sw*f_dot;
113 : }
114 :
115 : }
116 : }
117 :
|
