1 #include "ExternalMagneticFieldExtendedForce.h" 2 #include "Vector3DBlock.h" 3 #include "GenericTopology.h" 4 #include "Parallel.h" 5 #include "ScalarStructure.h" 6 7 using std::vector; 8 using std::string; 9 using namespace ProtoMol::Report; 10 11 namespace ProtoMol { 12 //_________________________________________________________________ ExternalMagneticFieldExtendedForce 13 14 15 16 const string ExternalMagneticFieldExtendedForce::keyword("ExternalMagneticField"); 17 static const Real SI_FORCE_FACTOR = Constant::SI::ELECTRON_CHARGE*Constant::SI::AVOGADRO*Constant::SI::KCAL/Constant::SI::LENGTH_AA*(1e8/sqrt(4184.0)); 18 ExternalMagneticFieldExtendedForce()19 ExternalMagneticFieldExtendedForce::ExternalMagneticFieldExtendedForce():ExtendedForce(),myB(Vector3D(0.0,0.0,0.0)){} 20 21 ExternalMagneticFieldExtendedForce(const Vector3D & b)22 ExternalMagneticFieldExtendedForce::ExternalMagneticFieldExtendedForce(const Vector3D& b ):ExtendedForce(),myB(b){} 23 24 getParameters(vector<Parameter> & parameters) const25 void ExternalMagneticFieldExtendedForce::getParameters(vector<Parameter>& parameters) const { 26 parameters.push_back(Parameter("-B",Value(myB),Text("External magnetic field [T]"))); 27 } 28 29 doMake(string &,vector<Value> values) const30 Force* ExternalMagneticFieldExtendedForce::doMake(string&, vector<Value> values) const { 31 return new ExternalMagneticFieldExtendedForce(values[0]); 32 } 33 evaluate(const GenericTopology * topo,const Vector3DBlock * positions,const Vector3DBlock * velocities,Vector3DBlock * forces,ScalarStructure * energies)34 void ExternalMagneticFieldExtendedForce::evaluate(const GenericTopology* topo, 35 const Vector3DBlock* positions, 36 const Vector3DBlock* velocities, 37 Vector3DBlock* forces, 38 ScalarStructure* energies){ 39 doEvaluate(topo,positions,velocities,forces,energies,0,positions->size()); 40 } 41 42 parallelEvaluate(const GenericTopology * topo,const Vector3DBlock * positions,const Vector3DBlock * velocities,Vector3DBlock * forces,ScalarStructure * energies)43 void ExternalMagneticFieldExtendedForce::parallelEvaluate(const GenericTopology* topo, 44 const Vector3DBlock* positions, 45 const Vector3DBlock* velocities, 46 Vector3DBlock* forces, 47 ScalarStructure* energies){ 48 unsigned int n = positions->size(); 49 unsigned int count = numberOfBlocks(topo,positions); 50 51 for(unsigned int i = 0;i<count;i++){ 52 if(Parallel::next()){ 53 int to = (n*(i+1))/count; 54 if(to > (int)n) 55 to = n; 56 int from = (n*i)/count; 57 doEvaluate(topo,positions,velocities,forces,energies,from,to); 58 } 59 } 60 } 61 62 doEvaluate(const GenericTopology * topo,const Vector3DBlock * positions,const Vector3DBlock * velocities,Vector3DBlock * forces,ScalarStructure *,int from,int to)63 void ExternalMagneticFieldExtendedForce::doEvaluate(const GenericTopology* topo, 64 const Vector3DBlock* positions, 65 const Vector3DBlock* velocities, 66 Vector3DBlock* forces, 67 ScalarStructure* /*energies*/, 68 int from, int to){ 69 Vector3D b(myB*SI_FORCE_FACTOR); 70 Real e = 0.0; 71 for(int i=from;i<to;i++){ 72 Vector3D f((*velocities)[i].cross(b)*topo->atomTypes[topo->atoms[i].type].charge); 73 //report << (*forces)[i] <<","; 74 (*forces)[i] += f; 75 //report << b <<","<< SI_FORCE_FACTOR<<","; 76 //report << (*velocities)[i].cross(b)*topo->atomTypes[topo->atoms[i].type].charge <<endr; 77 e += f*(*positions)[i]; 78 } 79 //(*energies)[ScalarStructure::OTHER] += e; 80 81 } 82 83 numberOfBlocks(const GenericTopology *,const Vector3DBlock * pos)84 unsigned int ExternalMagneticFieldExtendedForce::numberOfBlocks(const GenericTopology*, 85 const Vector3DBlock* pos){ 86 return std::min(static_cast<unsigned int>(Parallel::getAvailableNum()),static_cast<unsigned int>(pos->size())); 87 } 88 89 } 90