Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 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 12 : 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 12 : 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 12 : 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 12 : keys.add("numbered","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
45 : "The following provides information on the \\ref switchingfunction that are available. "
46 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
47 3 : }
48 :
49 1 : AlignedMatrixBase::AlignedMatrixBase( const ActionOptions& ao ):
50 : Action(ao),
51 1 : AdjacencyMatrixBase(ao)
52 : {
53 : // Read in the atomic positions
54 4 : readMaxTwoSpeciesMatrix( "ATOMS", "ATOMSA", "ATOMSB", true );
55 1 : unsigned nrows, ncols; retrieveTypeDimensions( nrows, ncols, ncol_t );
56 1 : if( mybasemulticolvars.size()==0 ) error("cannot use atom indices as input to this variable / input was not specified");
57 1 : if( getSizeOfInputVectors()<3 ) error("base multicolvars do not calculate an orientation");
58 : // Read in the switching function
59 1 : switchingFunction.resize( nrows, ncols );
60 2 : parseConnectionDescriptions("SWITCH",false,ncol_t);
61 :
62 : // Find the largest sf cutoff
63 1 : double sfmax=switchingFunction(0,0).get_dmax();
64 3 : for(unsigned i=0; i<getNumberOfNodeTypes(); ++i) {
65 3 : for(unsigned j=0; j<getNumberOfNodeTypes(); ++j) {
66 1 : double tsf=switchingFunction(i,j).get_dmax();
67 1 : if( tsf>sfmax ) sfmax=tsf;
68 : }
69 : }
70 : // And set the link cell cutoff
71 1 : setLinkCellCutoff( sfmax );
72 1 : }
73 :
74 2 : void AlignedMatrixBase::setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc ) {
75 2 : plumed_assert( id<2 );
76 2 : if( id==0 ) {
77 2 : plumed_assert( desc.size()==1 ); std::string errors; switchingFunction(j,i).set(desc[0],errors);
78 1 : if( errors.length()!=0 ) error("problem reading switching function in SWITCH keywrd description " + errors);
79 1 : if( j!=i) switchingFunction(i,j).set(desc[0],errors);
80 4 : log.printf(" %u th and %u th multicolvar groups must be within %s\n",i+1,j+1,(switchingFunction(i,j).description()).c_str() );
81 1 : } else if( id==1 ) {
82 1 : readOrientationConnector( i, j, desc );
83 : }
84 2 : }
85 :
86 57063 : double AlignedMatrixBase::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
87 114126 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
88 171189 : if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).get_dmax2() ) return 1.0;
89 51257 : return 0.0;
90 : }
91 :
92 5806 : double AlignedMatrixBase::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
93 5806 : unsigned ncomp=getSizeOfInputVectors(); Vector ddistance;
94 5806 : std::vector<double> orient0(ncomp), orient1(ncomp), dorient0(ncomp), dorient1(ncomp);
95 11612 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
96 5806 : getInputData( 0, true, myatoms, orient0 ); getInputData( 1, true, myatoms, orient1 );
97 17418 : double f_dot = computeVectorFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t,
98 11612 : distance, orient0, orient1, ddistance, dorient0, dorient1 );
99 :
100 : // Retrieve the weight of the connection
101 11612 : double dfunc, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) )-ncol_t ).calculate( distance.modulo(), dfunc );
102 :
103 5806 : if( !doNotCalculateDerivatives() ) {
104 0 : addAtomDerivatives( 1, 0, (-dfunc)*f_dot*distance - sw*ddistance, myatoms );
105 0 : addAtomDerivatives( 1, 1, (+dfunc)*f_dot*distance + sw*ddistance, myatoms );
106 0 : myatoms.addBoxDerivatives( 1, (-dfunc)*f_dot*Tensor(distance,distance) - sw*extProduct( distance, ddistance ) );
107 :
108 : // Add derivatives of orientation
109 0 : for(unsigned k=2; k<orient0.size(); ++k) { dorient0[k]*=sw; dorient1[k]*=sw; }
110 0 : mergeInputDerivatives( 1, 2, orient0.size(), 0, dorient0, getInputDerivatives( 0, true, myatoms ), myatoms );
111 0 : mergeInputDerivatives( 1, 2, orient1.size(), 1, dorient1, getInputDerivatives( 1, true, myatoms ), myatoms );
112 : }
113 11612 : return sw*f_dot;
114 : }
115 :
116 : }
117 4839 : }
118 :
|
