Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "VectorMultiColvar.h"
23 : #include "multicolvar/BridgedMultiColvarFunction.h"
24 :
25 : namespace PLMD {
26 : namespace crystallization {
27 :
28 32 : void VectorMultiColvar::registerKeywords( Keywords& keys ) {
29 32 : MultiColvarBase::registerKeywords( keys );
30 32 : }
31 :
32 26 : VectorMultiColvar::VectorMultiColvar(const ActionOptions& ao):
33 : Action(ao),
34 : MultiColvarBase(ao),
35 26 : store_director(true)
36 : {
37 : setLowMemOption(true);
38 26 : }
39 :
40 26 : void VectorMultiColvar::setVectorDimensionality( const unsigned& ncomp ) {
41 26 : ncomponents = ncomp; resizeFunctions(); // This resize needs to be here to ensure buffers are set to correct size in base
42 26 : }
43 :
44 0 : void VectorMultiColvar::doNotCalculateDirector() {
45 0 : store_director=false; // Need a sanity check in here so that you don't use the same instance of Q4 to calcualte vectors and directors
46 0 : }
47 :
48 100463 : double VectorMultiColvar::compute( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const {
49 : // Now calculate the vector
50 100463 : calculateVector( myatoms );
51 : // Sort out the active derivatives
52 100463 : updateActiveAtoms( myatoms );
53 :
54 : // Now calculate the norm of the vector (this is what we return here)
55 : double norm=0;
56 3888947 : for(unsigned i=0; i<ncomponents; ++i) norm += myatoms.getValue(2+i)*myatoms.getValue(2+i);
57 100463 : norm=sqrt(norm);
58 :
59 100463 : if( !doNotCalculateDerivatives() ) {
60 74901 : double inorm = 1.0 / norm; std::vector<double> dervec( ncomponents );
61 4344771 : for(unsigned i=0; i<ncomponents; ++i) dervec[i] = inorm*myatoms.getValue(2+i);
62 :
63 : MultiValue& myvals=myatoms.getUnderlyingMultiValue();
64 15502599 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
65 7713849 : unsigned jder=myvals.getActiveIndex(j);
66 533955495 : for(unsigned i=0; i<ncomponents; ++i) myvals.addDerivative( 1, jder, dervec[i]*myvals.getDerivative( 2+i, jder ) );
67 : }
68 : }
69 :
70 100463 : return norm;
71 : }
72 :
73 64156 : void VectorMultiColvar::normalizeVector( std::vector<double>& vals ) const {
74 : double inorm = 1.0;
75 64156 : if( vals[1]>epsilon ) inorm = 1.0 / vals[1];
76 4832132 : for(unsigned i=2; i<vals.size(); ++i) vals[i] = inorm*vals[i];
77 64156 : }
78 :
79 29195 : void VectorMultiColvar::normalizeVectorDerivatives( MultiValue& myvals ) const {
80 29195 : double v = myvals.get(1), weight = 1.0 / v, wdf = 1.0 / ( v*v*v );
81 7112123 : for(unsigned j=0; j<myvals.getNumberActive(); ++j) {
82 : double comp2=0.0; unsigned jder=myvals.getActiveIndex(j);
83 168107856 : for(unsigned jcomp=2; jcomp<myvals.getNumberOfValues(); ++jcomp) comp2 += myvals.get(jcomp)*myvals.getDerivative( jcomp, jder );
84 168107856 : for(unsigned jcomp=2; jcomp<myvals.getNumberOfValues(); ++jcomp) {
85 164566392 : myvals.setDerivative( jcomp, jder, weight*myvals.getDerivative( jcomp, jder ) - wdf*comp2*myvals.get(jcomp) );
86 : }
87 : }
88 29195 : }
89 :
90 0 : void VectorMultiColvar::addForcesOnAtoms( const std::vector<double>& inforces ) {
91 : plumed_dbg_assert( inforces.size()==getNumberOfDerivatives() );
92 0 : std::vector<double> oldforces( getNumberOfDerivatives() );
93 0 : getForcesFromVessels( oldforces );
94 0 : for(unsigned i=0; i<getNumberOfDerivatives(); ++i) oldforces[i]+=inforces[i];
95 0 : setForcesOnAtoms( oldforces );
96 0 : }
97 :
98 : }
99 4839 : }
|