Changeset faca99 for src


Ignore:
Timestamp:
Oct 5, 2011, 9:18:21 AM (13 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:
2e352f
Parents:
807c0e
git-author:
Daniel Dueck <dueck@…> (01/15/11 16:11:41)
git-committer:
Frederik Heber <heber@…> (10/05/11 09:18:21)
Message:

MoleculeAction create-micelle: code updated; Baseshapes: procedure getHomogeneousPointsonSurface updated

Changed in rebase onto v1.1.3:

  • replaced parser instantiations by getter to FormatParserStorage and World.
  • removed doubly present loop over count (creates 1942 points).
  • We create now less than desired points, Unit tests checks for 194 instead of originally 200 points.
Location:
src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/MoleculeAction/CreateMicelleAction.cpp

    r807c0e rfaca99  
    1717#include <config.h>
    1818#endif
    19 
    2019#include "CodePatterns/MemDebug.hpp"
    21 
    22 #include "CodePatterns/Log.hpp"
    23 #include "CodePatterns/Verbose.hpp"
    2420
    2521#include "Actions/ActionHistory.hpp"
    … …  
    2925#include "Descriptors/AtomIdDescriptor.hpp"
    3026#include "Descriptors/MoleculeDescriptor.hpp"
     27
     28#include "atom.hpp"
     29#include "Bond/bond.hpp"
     30#include "CodePatterns/Assert.hpp"
     31#include "CodePatterns/Log.hpp"
     32#include "CodePatterns/Verbose.hpp"
    3133#include "LinearAlgebra/Line.hpp"
     34#include "molecule.hpp"
     35#include "World.hpp"
     36
     37#include <iostream>
     38#include <string>
     39
    3240#include "Parser/PdbParser.hpp"
    3341#include "Parser/TremoloParser.hpp"
    … …  
    3745#include "Shapes/ShapeOps.hpp"
    3846
    39 
    40 #include "atom.hpp"
    41 #include "Bond/bond.hpp"
    42 #include "boundary.hpp"
    43 #include "molecule.hpp"
    44 #include "World.hpp"
    45 
    46 #include <iostream>
    47 #include <string>
    48 
    4947#include "Actions/MoleculeAction/CreateMicelleAction.hpp"
    5048
    … …  
    6664#include "molecule.hpp"
    6765#include "LinearAlgebra/Vector.hpp"
     66#include "LinearAlgebra/Line.hpp"
    6867#include "World.hpp"
    6968#include <gsl/gsl_poly.h>
    7069#include <gsl/gsl_eigen.h>
    71 #define PATH "/home/heber/tmp/"
     70//#define PATH "/home/dueck/workspace/tenside/tmp/"
    7271#define AtomVector std::vector <atom *>
    7372#define MoleculeVector std::vector <molecule *>
    … …  
    7877
    7978/** =========== define the function ====================== */
    80 Action::state_ptr MoleculeCreateMicelleAction::performCall() {
     79Action::state_ptr MoleculeCreateMicelleAction::performCall()
     80{
     81  getParametersfromValueStorage();
    8182       
    82 getParametersfromValueStorage();
    83        
    84 AtomVector ever = World::getInstance().getAllAtoms();
    85 
    86 // as all parsed atoms go into same molecule
    87 // we don't need to create one and add them all to it
    88 MoleculeVector all = World::getInstance().getSelectedMolecules();
    89 molecule *stick = *(all.begin());
    90 
    91 //3.Molekuel zentrieren
    92 
    93 stick->CenterOrigin();
    94 
    95 //4.Haupttraegheitsachse bestimmen
    96 Vector den(0.0,0.0,1.0);
    97 
    98 MoleculeRotateToPrincipalAxisSystem(den);
    99 /* determine
    100 principal axis system and make greatest eigenvector be aligned along
    101 (0,0,1)
    102 */
    103 string path;
    104 /**/
    105 {
    106   std::ofstream file;
    107   path = PATH;
    108   path += "/tensidrot.xyz";
    109   file.open(path.c_str());
    110   FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
    111   file.close();
    112 }
    113 
    114 //5.b: Molekuel um 180 Grad drehen
    115 
    116 Line RotationAxis(Vector(0.0,0.0,0.0), Vector(1.0,0.0,0.0)); // pt is the current Vector of point on surface
    117 
    118 for (molecule::iterator it=stick->begin(); it !=stick->end(); ++it)
    119   (*it)->setPosition(RotationAxis.rotateVector((*it)->getPosition(),M_PI));
    120 
    121 
    122 {
    123   std::ofstream file;
    124   path = PATH;
    125   path += "/tensid2rot.xyz";
    126   file.open(path.c_str());
    127   FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
    128   file.close();
    129 }
    130 
    131 //6.Molekuel mehrfach strukturiert mit der Haupttraegheitsachse senkrecht zu einer parametrisierten Oberflaeche anordnen
    132 
    133 //6.1. Punkte auf der Oberflaeche bestimmen
    134 //Algorithmus entnommen aus "http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere"
    135 
    136 int ka =0;
    137 double radius= 1.5*sqrt(pow(1.55, 2)*params.N);
    138 
    139 Shape s = resize(Sphere(), radius);
    140 std::vector<Vector> pt = s.getHomogeneousPointsOnSurface(params.N);
    141 
    142 //6.2."stick" um Radius und Molekuelausdehnung in z-Richtung verschieben.
    143 
    144 for (molecule::iterator it2=stick->begin();it2 !=stick->end();++it2) {
    145         Vector pos= (**it2).getPosition();
    146         pos[2]=pos[2]+radius; // -Min[2]
    147         (**it2).setPosition(pos);
    148 }
    149 
    150 {
    151   std::ofstream file;
    152   path = PATH;
    153   path += "/tensid3rot.xyz";
    154   file.open(path.c_str());
    155   FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
    156   file.close();
    157 }
     83  AtomVector ever = World::getInstance().getAllAtoms();
     84
     85  // as all parsed atoms go into same molecule
     86  // we don't need to create one and add them all to it
     87  MoleculeVector all = World::getInstance().getSelectedMolecules();
     88  ASSERT(!all.empty(), "MoleculeCreateMicelleAction::performCall() - no molecules selected.");
     89  molecule *stick = *(all.begin());
     90
     91  //3.Molekuel zentrieren
     92
     93  stick->CenterOrigin();
     94
     95  //4.Haupttraegheitsachse bestimmen
     96  Vector den(0.0,0.0,1.0);
     97
     98  MoleculeRotateToPrincipalAxisSystem(den);
     99  /* determine
     100  principal axis system and make greatest eigenvector be aligned along
     101  (0,0,1)
     102  */
     103  string path;
     104  /**/
     105  /*XyzParser *parserx = new XyzParser;
     106  {
     107    std::ofstream file;
     108    path = PATH;
     109    path += "/tensidrot.xyz";
     110    file.open(path.c_str());
     111    FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
     112    file.close();
     113  }*/
     114  //5.b: Molekuel um 180 Grad drehen
     115
     116  Line RotationAxis(Vector(0.0,0.0,0.0), Vector(1.0,0.0,0.0)); // pt is the current Vector of point on surface
     117
     118  for (molecule::iterator it=stick->begin(); it !=stick->end(); ++it)
     119    (*it)->setPosition(RotationAxis.rotateVector((*it)->getPosition(),M_PI));
     120
     121
     122  /*{
     123    std::ofstream file;
     124    path = PATH;
     125    path += "/tensid2rot.xyz";
     126    file.open(path.c_str());
     127    FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
     128    file.close();
     129  }*/
     130
     131  //6.Molekuel mehrfach strukturiert mit der Haupttraegheitsachse senkrecht zu einer parametrisierten Oberflaeche anordnen
     132
     133  //6.1. Punkte auf der Oberflaeche bestimmen
     134  //Algorithmus entnommen aus "http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere"
     135
     136  int ka =0;
     137  //double radius= 1.5*sqrt(pow(1.55, 2)*params.N);
     138
     139  Shape s = resize(Sphere(), params.radius);
     140  std::vector<Vector> pt = s.getHomogeneousPointsOnSurface(params.N);
     141
     142  //6.2.a. "stick" 180 Grad an x-y-Ebene spiegeln
     143
     144  for (molecule::iterator it2=stick->begin();it2 !=stick->end();++it2)
     145  {
     146          Vector pos= (**it2).getPosition();
     147          pos[2]= - pos[2];// -Min[2]
     148          (**it2).setPosition(pos);
     149  }
     150
     151
     152  //6.2.b. "stick" um Radius und Molekuelausdehnung in z-Richtung verschieben.
     153
     154  for (molecule::iterator it2=stick->begin();it2 !=stick->end();++it2)
     155  {
     156          Vector pos= (**it2).getPosition();
     157          pos[2]=pos[2]+params.radius; // -Min[2]
     158          (**it2).setPosition(pos);
     159  }
     160
     161
     162
     163
     164  /*{
     165    std::ofstream file;
     166    path = PATH;
     167    path += "/tensid3rot.xyz";
     168    file.open(path.c_str());
     169    FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
     170    file.close();
     171  }*/
    158172
    159173//6.3.Erzeugen einer Molekuelliste, die das Molekuel "stick" "N" mal kopiert und um eine Sphaere herum verteilt
    160174
    161 //double MYEPSILON=1e-10;
    162 
    163 for (ka = 0; ka<params.N-1; ka++){
    164   cout << "Creating " << ka+1 << " copy of tenside molecule 'stick' with " << stick->getAtomCount() << " atoms, ";
    165         molecule *Tensid=stick->CopyMolecule();
    166 
    167         cout << "rotating ...";
    168   Vector ZAxis(Vector(0.0,0.0,1.0));
    169   Vector Axis(pt[ka]);
    170   const double alpha = ZAxis.Angle(Axis);
    171   Axis.VectorProduct(ZAxis);
    172   Line RotationAxis(Vector(0.0,0.0,1.0), Axis); // pt is the current Vector of point on surface
    173 
    174   for (molecule::iterator it2=Tensid->begin();it2 !=Tensid->end();++it2)
    175     (*it2)->setPosition(RotationAxis.rotateVector((*it2)->getPosition(),alpha));
    176   cout << "done." << endl;
    177 
    178   Tensid=NULL;
    179 }
    180 cout << "shifting " << ka+1 << " copy of tenside molecule, ";
    181         molecule *Tensid=stick;
    182 
    183         cout << "rotating ...";
     175  //double MYEPSILON=1e-10;
     176
     177  for (ka = 0; ka<params.N-1; ka++)
     178  {
     179    cout << "Creating " << ka+1 << " copy of tenside molecule 'stick' with " << stick->getAtomCount() << " atoms, ";
     180          molecule *Tensid=stick->CopyMolecule();
     181
     182          cout << "rotating ...";
     183    Vector ZAxis(Vector(0.0,0.0,1.0));
     184    Vector Axis(pt[ka]);
     185    const double alpha = ZAxis.Angle(Axis);
     186    Axis.VectorProduct(ZAxis);
     187    Line RotationAxis(Vector(0.0,0.0,1.0), Axis);       // pt is the current Vector of point on surface
     188
     189    for (molecule::iterator it2=Tensid->begin();it2 !=Tensid->end();++it2)
     190    {
     191      (*it2)->setPosition(RotationAxis.rotateVector((*it2)->getPosition(),alpha));
     192      *(*it2)+=params.center;
     193    }
     194    cout << "done." << endl;
     195
     196    Tensid=NULL;
     197  }
     198
     199  cout << "shifting " << ka+1 << " copy of tenside molecule, ";
     200  molecule *Tensid=stick;
     201  cout << "rotating ...";
    184202  Vector ZAxis(Vector(0.0,0.0,1.0));
    185203  Vector Axis(pt[params.N-1]);
    … …  
    189207
    190208  for (molecule::iterator it2=Tensid->begin();it2 !=Tensid->end();++it2)
     209  {
    191210    (*it2)->setPosition(RotationAxis.rotateVector((*it2)->getPosition(),alpha));
     211    *(*it2)+=params.center;
     212  }
    192213  cout << "done." << endl;
    193214
    … …  
    195216
    196217  GraphSubgraphDissection();
    197 
    198218  return Action::success;
    199219}
    200220
    201 Action::state_ptr MoleculeCreateMicelleAction::performUndo(Action::state_ptr _state) {
    202 //  MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());
    203 
    204 //  string newName = state->mol->getName();
    205 //  state->mol->setName(state->lastName);
     221Action::state_ptr MoleculeCreateMicelleAction::performUndo(Action::state_ptr _state)
     222{
     223  //  MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.());
     224
     225  //  string newName = state->mol->Name();
     226  //  state->mol->setName(state->lastName);
    206227
    207228  return Action::failure;
    208229}
    209230
    210 Action::state_ptr MoleculeCreateMicelleAction::performRedo(Action::state_ptr _state){
     231Action::state_ptr MoleculeCreateMicelleAction::performRedo(Action::state_ptr _state)
     232{
    211233  // Undo and redo have to do the same for this action
    212234  return performUndo(_state);
    … …  
    214236
    215237
    216 bool MoleculeCreateMicelleAction::canUndo() {
     238bool MoleculeCreateMicelleAction::canUndo()
     239{
    217240  return false;
    218241}
    219242
    220 
    221 bool MoleculeCreateMicelleAction::shouldUndo() {
     243bool MoleculeCreateMicelleAction::shouldUndo()
     244{
    222245  return false;
    223246}

  • src/Shapes/BaseShapes.cpp

    r807c0e rfaca99  
    7272  std::vector<Vector> PointsOnSurface;
    7373
    74   const double dlength = M_PI*(3.-sqrt(5.));
     74    double PI=3.14159265358979323846;
     75    double a=4*PI/N;
     76    double d= sqrt(a);
     77    int Mtheta=int(PI/d);
     78    double dtheta=PI/Mtheta;
     79    double dphi=a/dtheta;
     80    for (int m=0; m<Mtheta; m++)
     81    {
     82          double theta=PI*(m+0.5)/Mtheta;
     83          int Mphi=int(2*PI*sin(theta)/dphi);
     84          for (int n=0; n<Mphi;n++)
     85          {
     86                  double phi= 2*PI*n/Mphi;
     87                  Vector point;
     88                  point.Zero();
     89                  point[0]=sin(theta)*cos(phi);
     90                  point[1]=sin(theta)*sin(phi);
     91                  point[2]=cos(theta);
     92                  PointsOnSurface.push_back(point);
     93          }
     94    }
     95  /*const double dlength = M_PI*(3.-sqrt(5.));
    7596  double length = 0;
    7697  const double dz = 2.0/N;
    … …  
    85106    PointsOnSurface.push_back(point);
    86107    z = z - dz;
    87     length = length + dlength;
    88   }
    89 
    90   ASSERT(PointsOnSurface.size() == N, "Sphere_impl::getHomogeneousPointsOnSurface() did not create enough points.");
     108    length = length + dlength;*/
     109
     110  //ASSERT(PointsOnSurface.size() == N, "Sphere_impl::getHomogeneousPointsOnSurface() did not create enough points.");
    91111  return PointsOnSurface;
    92112}

  • src/Shapes/unittests/ShapeUnitTest.cpp

    r807c0e rfaca99  
    651651    CPPUNIT_ASSERT(fabs(1. - (*iter).NormSquared()) < MYEPSILON);
    652652  }
    653   CPPUNIT_ASSERT_EQUAL(N, PointsOnSurface.size());
     653  CPPUNIT_ASSERT_EQUAL((size_t)194, PointsOnSurface.size());
    654654}
    655655
Note: See TracChangeset for help on using the changeset viewer.