Changeset 4e855e for src/molecule_dynamics.cpp
- Timestamp:
- Feb 4, 2011, 6:56:38 PM (15 years ago)
- Branches:
- Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
- Children:
- febef3
- Parents:
- fd3788
- git-author:
- Frederik Heber <heber@…> (02/04/11 18:54:12)
- git-committer:
- Frederik Heber <heber@…> (02/04/11 18:56:38)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/molecule_dynamics.cpp
rfd3788 r4e855e 574 574 { 575 575 Info FunctionInfo(__func__); 576 ifstream input(file);577 576 string token; 578 577 stringstream item; … … 582 581 583 582 const int AtomCount = getAtomCount(); 584 // check file 585 if (input == NULL) { 583 // parse file into ForceMatrix 584 std::ifstream input(file); 585 if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) { 586 DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl); 587 performCriticalExit(); 586 588 return false; 587 } else { 588 // parse file into ForceMatrix 589 if (!Force.ParseMatrix(file, 0,0,0)) { 590 DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl); 591 performCriticalExit(); 592 return false; 593 } 594 if (Force.RowCounter[0] != AtomCount) { 595 DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl); 596 performCriticalExit(); 597 return false; 598 } 599 // correct Forces 600 Velocity.Zero(); 601 for(int i=0;i<AtomCount;i++) 602 for(int d=0;d<NDIM;d++) { 603 Velocity[d] += Force.Matrix[0][i][d+offset]; 604 } 605 for(int i=0;i<AtomCount;i++) 606 for(int d=0;d<NDIM;d++) { 607 Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount); 608 } 609 // solve a constrained potential if we are meant to 610 if (configuration.DoConstrainedMD) { 611 // calculate forces and potential 612 atom **PermutationMap = NULL; 613 ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem()); 614 EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force); 615 delete[](PermutationMap); 616 } 617 618 // and perform Verlet integration for each atom with position, velocity and force vector 619 // check size of vectors 620 for_each(atoms.begin(), 621 atoms.end(), 622 boost::bind(&atom::VelocityVerletUpdate,_1,MDSteps+1, &configuration, &Force, (const size_t) 0)); 623 } 589 } 590 input.close(); 591 if (Force.RowCounter[0] != AtomCount) { 592 DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl); 593 performCriticalExit(); 594 return false; 595 } 596 // correct Forces 597 Velocity.Zero(); 598 for(int i=0;i<AtomCount;i++) 599 for(int d=0;d<NDIM;d++) { 600 Velocity[d] += Force.Matrix[0][i][d+offset]; 601 } 602 for(int i=0;i<AtomCount;i++) 603 for(int d=0;d<NDIM;d++) { 604 Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount); 605 } 606 // solve a constrained potential if we are meant to 607 if (configuration.DoConstrainedMD) { 608 // calculate forces and potential 609 atom **PermutationMap = NULL; 610 ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem()); 611 EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force); 612 delete[](PermutationMap); 613 } 614 615 // and perform Verlet integration for each atom with position, velocity and force vector 616 // check size of vectors 617 for_each(atoms.begin(), 618 atoms.end(), 619 boost::bind(&atom::VelocityVerletUpdate,_1,MDSteps+1, &configuration, &Force, (const size_t) 0)); 620 624 621 // correct velocities (rather momenta) so that center of mass remains motionless 625 622 Velocity = atoms.totalMomentumAtStep(MDSteps+1);
Note:
See TracChangeset
for help on using the changeset viewer.