Changeset 2440ce

Timestamp:

Jan 14, 2015, 9:03:24 PM (10 years ago)

Author:

Frederik Heber <heber@…>

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:

cad383

Parents:

601ef8

git-author:

Frederik Heber <heber@…> (01/14/15 20:31:24)

git-committer:

Frederik Heber <heber@…> (01/14/15 21:03:24)

Message:

FIX: Changed SuspendInMoleculeAction to catch segfault when rho=1 was given.

• however, the action is still not tested to work.
• TESTFIX: suspend-in-water gets parameter density via "density" not via "suspend-in-water", set desired density to faulting 1 for the moment.
• TESTFIX: added water molecule to test.conf, to have at least two molecules.

Files:

• src/Actions/FillAction/SuspendInMoleculeAction.cpp (modified) (6 diffs)
• src/Actions/FillAction/SuspendInMoleculeAction.def (modified) (2 diffs)
• tests/regression/Filling/SuspendInWater/pre/single_molecule.conf (added)
• tests/regression/Filling/SuspendInWater/pre/test.conf (modified) (5 diffs)
• tests/regression/Filling/SuspendInWater/testsuite-suspend-in-water.at (modified) (4 diffs)

src/Actions/FillAction/SuspendInMoleculeAction.cpp

@@ -129 +129 @@
   filler->CenterEdge();
 
+  std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
+  if (molecules.size() < 2) {
+    STATUS("There must be at least two molecules: filler and to be suspended.");
+    return Action::failure;
+  }
+
   /// first we need to calculate some volumes and masses
-  std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
   double totalmass = 0.;
   const bool IsAngstroem = true;
@@ -139 +144 @@
       iter != molecules.end(); ++iter)
   {
+    // skip the filler
+    if (*iter == filler)
+      continue;
     molecule & mol = **iter;
     const double mass = calculateMass(mol);
@@ -152 +160 @@
   LOG(1, "INFO: The average density is " << setprecision(10)
       << totalmass / clustervolume << " atomicmassunit/"
-      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
+      << (IsAngstroem ? " angstrom" : " atomiclength") << "^3.");
+  if ( ((totalmass / clustervolume < 1.) && (params.density.get() > 1.))
+      || ((totalmass / clustervolume > 1.) && (params.density.get() < 1.))) {
+    STATUS("Desired and present molecular densities must both be either in [0,1) or in (1, inf).");
+    return Action::failure;
+  }
 
   // calculate maximum solvent density
   std::vector<double> fillerdiameters(NDIM, 0.);
-  const double solventdensity =
-      calculateMass(*filler) / calculateEnvelopeVolume(*filler, fillerdiameters);
-
-  std::vector<unsigned int> counts(3, 0);
-  Vector offset(.5,.5,.5);
+  const double fillervolume = calculateEnvelopeVolume(*filler, fillerdiameters);
+  const double fillermass = calculateMass(*filler);
+  LOG(1, "INFO: The filler's mass is " << setprecision(10)
+      << fillermass << " atomicmassunit, and it's volume is "
+      << fillervolume << (IsAngstroem ? " angstrom" : " atomiclength") << "^3.");
+  const double solventdensity = fillermass / fillervolume;
 
   /// solve cubic polynomial
   double cellvolume = 0.;
   LOG(1, "Solving equidistant suspension in water problem ...");
-  cellvolume = (totalmass / solventdensity
-      - (totalmass / clustervolume)) / (params.density.get() - 1.);
+  // s = solvent, f = filler, 0 = initial molecules/cluster
+  // v_s = v_0 + v_f, m_s = m_0 + rho_f * v_f --> rho_s = m_s/v_s ==> v_f = (m_0 - rho_s * v_o) / (rho_s - rho_f)
+  cellvolume = (totalmass - params.density.get() * clustervolume) / (params.density.get() - 1.) + clustervolume;
   LOG(1, "Cellvolume needed for a density of " << params.density.get()
       << " g/cm^3 is " << cellvolume << " angstroem^3.");
@@ -173 +188 @@
       (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
   LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is "
-      << minimumvolume << "angstrom^3.");
+      << minimumvolume << " angstrom^3.");
 
   if (minimumvolume > cellvolume) {
@@ -187 +202 @@
         + GreatestDiameter[1] * GreatestDiameter[2];
     BoxLengths[2] = minimumvolume - cellvolume;
-    double x0 = 0.;
-    double x1 = 0.;
-    double x2 = 0.;
+    std::vector<double> x(3, 0.);
     // for cubic polynomial there are either 1 or 3 unique solutions
-    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1)
-      LOG(0, "RESULT: The resulting spacing is: " << x0 << " .");
-    else {
-      LOG(0, "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " .");
-      std::swap(x0, x2); // sorted in ascending order
-    }
+    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x[0], &x[1], &x[2]) == 1) {
+      x[1] = x[0];
+      x[2] = x[0];
+    } else {
+      std::swap(x[0], x[2]); // sorted in ascending order
+    }
+    LOG(0, "RESULT: The resulting spacing is: " << x << " .");
 
     cellvolume = 1.;
     for (size_t i = 0; i < NDIM; ++i) {
-      BoxLengths[i] = x0 + GreatestDiameter[i];
+      BoxLengths[i] = x[i] + GreatestDiameter[i];
       cellvolume *= BoxLengths[i];
     }
-
-    // set new box dimensions
-    LOG(0, "Translating to box with these boundaries.");
-    {
-      RealSpaceMatrix domain;
-      for(size_t i =0; i<NDIM;++i)
-        domain.at(i,i) = BoxLengths[i];
-      World::getInstance().setDomain(domain);
-    }
-//    mol->CenterInBox();
+  }
+
+  // TODO: Determine counts from resulting mass correctly (hard problem due to integers)
+  std::vector<unsigned int> counts(3, 0);
+  const unsigned int totalcounts = round(params.density.get() * cellvolume - totalmass) / fillermass;
+  if (totalcounts > 0) {
+    counts[0] = ceil(BoxLengths[0]/3.1);
+    counts[1] = ceil(BoxLengths[1]/3.1);
+    counts[2] = ceil(BoxLengths[2]/3.1);
   }
 
@@ -231 +244 @@
       params.RandMoleculeDisplacement.get(),
       params.DoRotate.get());
+  Vector offset(.5,.5,.5);
   filler_preparator.addCubeMesh(
       counts,

src/Actions/FillAction/SuspendInMoleculeAction.def

@@ -17 +17 @@
 #include "Parameters/Validators/RangeValidator.hpp"
 #include "Parameters/Validators/STLVectorValidator.hpp"
+#include "Parameters/Validators/Ops_Validator.hpp"
 #include "Parameters/Validators/Specific/BoxLengthValidator.hpp"
 #include "Parameters/Validators/Specific/VectorZeroOneComponentsValidator.hpp"
@@ -25 +26 @@
 #define paramtypes (double)(double)(double)(double)(bool)
 #define paramtokens ("density")("min-distance")("random-atom-displacement")("random-molecule-displacement")("DoRotate")
-#define paramdescriptions ("desired density for the total domain")("minimum distance of water molecules to present atoms")("magnitude of random atom displacement")("magnitude of random molecule displacement")("whether to rotate or not")
+#define paramdescriptions ("desired density for the total domain, unequal 1.")("minimum distance of water molecules to present atoms")("magnitude of random atom displacement")("magnitude of random molecule displacement")("whether to rotate or not")
 #define paramdefaults (PARAM_DEFAULT(1.))(PARAM_DEFAULT(1.))(PARAM_DEFAULT(0.))(PARAM_DEFAULT(0.))(PARAM_DEFAULT(false))
 #define paramreferences (density)(mindistance)(RandAtomDisplacement)(RandMoleculeDisplacement)(DoRotate)
 #define paramvalids \
-(RangeValidator< double >(0., std::numeric_limits<double>::max())) \
+(RangeValidator< double >(0., 1. - std::numeric_limits<double>::epsilon()) || RangeValidator< double >(1. + std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::max())) \
 (BoxLengthValidator()) \
 (BoxLengthValidator()) \

tests/regression/Filling/SuspendInWater/pre/test.conf

@@ -28 +28 @@
 OutSrcStep      5       # Output "restart" data every ..th step
 TargetTemp      0.000950045     # Target temperature
-MaxPsiStep      0       # number of Minimisation steps per state (0 - default)
+MaxPsiStep      3       # number of Minimisation steps per state (0 - default)
 EpsWannier      1e-07   # tolerance value for spread minimisation of orbitals
 
@@ -35 +35 @@
 RelEpsTotalE    1e-07   # relative change in total energy
 RelEpsKineticE  1e-05   # relative change in kinetic energy
-MaxMinStopStep  0       # check every ..th steps
-MaxMinGapStopStep       0       # check every ..th steps
+MaxMinStopStep  12      # check every ..th steps
+MaxMinGapStopStep       1       # check every ..th steps
 
 # Values specifying when to stop for INIT, otherwise same as above
@@ -42 +42 @@
 InitRelEpsTotalE        1e-05   # relative change in total energy
 InitRelEpsKineticE      0.0001  # relative change in kinetic energy
-InitMaxMinStopStep      0       # check every ..th steps
-InitMaxMinGapStopStep   0       # check every ..th steps
+InitMaxMinStopStep      12      # check every ..th steps
+InitMaxMinGapStopStep   1       # check every ..th steps
 
 BoxLength                       # (Length of a unit cell)
@@ -55 +55 @@
 RiemannTensor   0       # (Use metric)
 PsiType         0       # 0 - doubly occupied, 1 - SpinUp,SpinDown
-MaxPsiDouble    0       # here: specifying both maximum number of SpinUp- and -Down-states
-PsiMaxNoUp      0       # here: specifying maximum number of SpinUp-states
-PsiMaxNoDown    0       # here: specifying maximum number of SpinDown-states
+MaxPsiDouble    12      # here: specifying both maximum number of SpinUp- and -Down-states
+PsiMaxNoUp      12      # here: specifying maximum number of SpinUp-states
+PsiMaxNoDown    12      # here: specifying maximum number of SpinDown-states
 AddPsis         0       # Additional unoccupied Psis for bandgap determination
 
@@ -64 +64 @@
 IsAngstroem     1       # 0 - Bohr, 1 - Angstroem
 RelativeCoord   0       # whether ion coordinates are relative (1) or absolute (0)
-MaxTypes        2       # maximum number of different ion types
+MaxTypes        3       # maximum number of different ion types
 
 # Ion type data (PP = PseudoPotential, Z = atomic number)
 #Ion_TypeNr.    Amount  Z       RGauss  L_Max(PP)L_Loc(PP)IonMass       # chemical name, symbol
-Ion_Type1       8       1       1.0     3       3       1.00800000000   Hydrogen        H
+Ion_Type1       10      1       1.0     3       3       1.00800000000   Hydrogen        H
 Ion_Type2       3       6       1.0     3       3       12.01100000000  Carbon  C
+Ion_Type3       1       8       1.0     3       3       15.99900000000  Oxygen  O
 #Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)
-Ion_Type2_1     9.782085945     3.275186040     3.535886037     0 # molecule nr 0
-Ion_Type2_2     8.532785963     4.158586027     3.535886037     0 # molecule nr 1
-Ion_Type2_3     7.283585982     3.275186040     3.535886037     0 # molecule nr 2
-Ion_Type1_1     9.782085945     2.645886050     2.645886050     0 # molecule nr 3
-Ion_Type1_2     9.782085945     2.645886050     4.425886024     0 # molecule nr 4
-Ion_Type1_3     10.672039608    3.904536878     3.535886037     0 # molecule nr 5
-Ion_Type1_4     8.532785963     4.787886018     2.645886050     0 # molecule nr 6
-Ion_Type1_5     8.532785963     4.787886018     4.425886024     0 # molecule nr 7
-Ion_Type1_6     6.393632318     3.904536877     3.535886037     0 # molecule nr 8
-Ion_Type1_7     7.283585982     2.645886050     2.645886050     0 # molecule nr 9
-Ion_Type1_8     7.283585982     2.645886050     4.425886024     0 # molecule nr 10
+Ion_Type1_1     9.782085945     2.645886050     2.645886050     0 # molecule nr 0
+Ion_Type1_2     9.782085945     2.645886050     4.425886024     0 # molecule nr 1
+Ion_Type1_3     10.672039608    3.904536878     3.535886037     0 # molecule nr 2
+Ion_Type1_4     8.532785963     4.787886018     2.645886050     0 # molecule nr 3
+Ion_Type1_5     8.532785963     4.787886018     4.425886024     0 # molecule nr 4
+Ion_Type1_6     6.393632318     3.904536877     3.535886037     0 # molecule nr 5
+Ion_Type1_7     7.283585982     2.645886050     2.645886050     0 # molecule nr 6
+Ion_Type1_8     7.283585982     2.645886050     4.425886024     0 # molecule nr 7
+Ion_Type2_1     9.782085945     3.275186040     3.535886037     0 # molecule nr 8
+Ion_Type2_2     8.532785963     4.158586027     3.535886037     0 # molecule nr 9
+Ion_Type2_3     7.283585982     3.275186040     3.535886037     0 # molecule nr 10
+Ion_Type3_1     2.000000000     2.000000000     2.000000000     0 # molecule nr 11
+Ion_Type1_9     2.758602000     2.000000000     2.504284000     0 # molecule nr 12
+Ion_Type1_10    2.758602000     2.000000000     1.495716000     0 # molecule nr 13

tests/regression/Filling/SuspendInWater/testsuite-suspend-in-water.at

@@ -18 +18 @@
 ### suspend in water with certain density
 
+AT_SETUP([Filling - suspend in water fails with rho=1])
+AT_KEYWORDS([filling suspend-in-water])
+
+file=test.conf
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/test.conf $file], 0)
+AT_CHECK([chmod u+w $file], 0)
+AT_CHECK([../../molecuilder -i $file -I -v 3 --select-molecule-by-id 0 -u --density 1.], 5, [stdout], [stderr])
+
+AT_CLEANUP
+
+AT_SETUP([Filling - suspend in water fails with just one molecule])
+AT_KEYWORDS([filling suspend-in-water])
+
+file=single_molecule.conf
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/$file $file], 0)
+AT_CHECK([chmod u+w $file], 0)
+AT_CHECK([../../molecuilder -i $file -I -v 3 --select-molecule-by-id 0 -u --density 1.], 5, [stdout], [stderr])
+AT_CHECK([grep "at least two molecules" stdout], 0, [ignore], [ignore])
+
+AT_CLEANUP
+
+# rho must be on same side of 1 as the present density
+AT_SETUP([Filling - suspend in water fails with wrong rho])
+AT_KEYWORDS([filling suspend-in-water])
+
+file=test.conf
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/$file $file], 0)
+AT_CHECK([chmod u+w $file], 0)
+AT_CHECK([../../molecuilder -i $file -I -v 3 --select-molecule-by-id 0 -u --density .5], 5, [stdout], [stderr])
+AT_CHECK([grep "Desired and present molecular densities must both be either" stdout], 0, [ignore], [ignore])
+
+AT_CLEANUP
+
 AT_SETUP([Filling - suspend in water])
 AT_XFAIL_IF([/bin/true])
@@ -23 +56 @@
 
 file=test.conf
-AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/test.conf $file], 0)
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/$file $file], 0)
 AT_CHECK([chmod u+w $file], 0)
-AT_CHECK([../../molecuilder -i $file  -v 3 --select-molecule-by-id 0 -u 1.3], 0, [stdout], [stderr])
+AT_CHECK([../../molecuilder -i $file -I -v 3 --select-molecule-by-id 0 -u --density 2.], 5, [stdout], [stderr])
 AT_CHECK([diff $file ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/post/test.conf], 0, [ignore], [ignore])
 
 AT_CLEANUP
-
 
 AT_SETUP([Filling - suspend in water with Undo])
@@ -38 +70 @@
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/test.conf $file], 0)
 AT_CHECK([chmod u+w $file], 0)
-AT_CHECK([../../molecuilder -i $file  -v 3 --select-molecule-by-id 0 -u 1.3 --undo], 0, [stdout], [stderr])
+AT_CHECK([../../molecuilder -i $file -I -v 3 --select-molecule-by-id 0 -u --density 2. --undo], 0, [stdout], [stderr])
 AT_CHECK([diff $file ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/post/test.conf], 0, [ignore], [ignore])
 
@@ -51 +83 @@
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/pre/test.conf $file], 0)
 AT_CHECK([chmod u+w $file], 0)
-AT_CHECK([../../molecuilder -i $file  -v 3 --select-molecule-by-id 0 -u 1.3 --undo --redo], 0, [stdout], [stderr])
+AT_CHECK([../../molecuilder -i $file -I -v 3 --select-molecule-by-id 0 -u --density 2. --undo --redo], 0, [stdout], [stderr])
 AT_CHECK([diff $file ${abs_top_srcdir}/tests/regression/Filling/SuspendInWater/post/test.conf], 0, [ignore], [ignore])