Changeset f7f7a4


Ignore:
Timestamp:
Sep 28, 2009, 6:41:52 PM (15 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:
d52e70c
Parents:
8cd903
git-author:
Frederik Heber <heber@…> (09/28/09 18:39:28)
git-committer:
Frederik Heber <heber@…> (09/28/09 18:41:52)
Message:

Implemenation of embedding merge, untested. LinearInterpolation now has switch for using identity map.

  • embed merge in MoleculeListClass determines the non-convex envelope and tries to fill the domain with the filler molecule everywhere outside of this domain
  • is added to MergeMolecule in builder.cpp, case 'e'
  • in molecule::LinearInterpolationBetweenConfiguration() now has additional parameter to tell whether we have to look for a mapping or use the identity.
Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • src/builder.cpp

    r8cd903 rf7f7a4  
    10171017      cin >> nr;
    10181018      count = 1;
    1019       for( MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1019      for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    10201020        if (nr == (*ListRunner)->IndexNr) {
    10211021          mol = *ListRunner;
    10221022          molecules->ListOfMolecules.erase(ListRunner);
    10231023          delete(mol);
     1024          break;
    10241025        }
    10251026      break;
     
    10741075
    10751076    case 'e':
    1076       cout << Verbose(0) << "Not implemented yet." << endl;
     1077      {
     1078        int src, dest;
     1079        molecule *srcmol = NULL, *destmol = NULL;
     1080        do {
     1081          cout << Verbose(0) << "Enter index of matrix molecule (the variable one): ";
     1082          cin >> src;
     1083          srcmol = molecules->ReturnIndex(src);
     1084        } while ((srcmol == NULL) && (src != -1));
     1085        do {
     1086          cout << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ";
     1087          cin >> dest;
     1088          destmol = molecules->ReturnIndex(dest);
     1089        } while ((destmol == NULL) && (dest != -1));
     1090        if ((src != -1) && (dest != -1))
     1091          molecules->EmbedMerge(destmol, srcmol);
     1092      }
    10771093      break;
    10781094
     
    16251641                cout << Verbose(0) << "Evaluating non-convex envelope.";
    16261642                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
    1627                 start = clock();
     1643                start = clock();
    16281644                LinkedCell LCList(mol, atof(argv[argptr])*2.);
    1629                 FindNonConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr+1], atof(argv[argptr]));
     1645                FindNonConvexBorder((ofstream *)&cout, mol, &LCList, atof(argv[argptr]), argv[argptr+1]);
    16301646                //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
    1631                 end = clock();
    1632                 cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     1647                end = clock();
     1648                cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    16331649                argptr+=2;
    16341650              }
     
    16531669            case 'L':
    16541670              ExitFlag = 1;
    1655               SaveFlag = true;
    1656               cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
    1657               if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration))
    1658                 cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
    1659               else
    1660                 cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
    1661               argptr+=3;
     1671              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1672                ExitFlag = 255;
     1673                cerr << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl;
     1674              } else {
     1675                SaveFlag = true;
     1676                cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
     1677                if (atoi(argv[argptr+3]) == 1)
     1678                  cout << Verbose(1) << "Using Identity for the permutation map." << endl;
     1679                if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false)
     1680                  cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
     1681                else
     1682                  cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
     1683                argptr+=4;
     1684              }
    16621685              break;
    16631686            case 'P':
     
    17091732                ExitFlag = 1;
    17101733                SaveFlag = true;
    1711                 cout << Verbose(1) << "Translating all ions to new origin." << endl;
     1734                cout << Verbose(1) << "Translating all ions by given vector." << endl;
    17121735                for (int i=NDIM;i--;)
    17131736                  x.x[i] = atof(argv[argptr+i]);
     
    17151738                argptr+=3;
    17161739              }
     1740              break;
    17171741            case 'T':
    17181742              ExitFlag = 1;
     
    17231747                ExitFlag = 1;
    17241748                SaveFlag = true;
    1725                 cout << Verbose(1) << "Translating all ions periodically to new origin." << endl;
     1749                cout << Verbose(1) << "Translating all ions periodically by given vector." << endl;
    17261750                for (int i=NDIM;i--;)
    17271751                  x.x[i] = atof(argv[argptr+i]);
     
    18641888            case 'o':
    18651889              ExitFlag = 1;
    1866               if ((argptr >= argc) || (argv[argptr][0] == '-')){
     1890              if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){
    18671891                ExitFlag = 255;
    1868                 cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
     1892                cerr << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;
    18691893              } else {
    18701894                cout << Verbose(0) << "Evaluating volume of the convex envelope.";
    1871                 cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
     1895                cout << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;
     1896                cout << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;
    18721897                LinkedCell LCList(mol, 10.);
    18731898                //FindConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr]);
    1874                 FindNonConvexBorder((ofstream *)&cout, mol, &LCList, argv[argptr], 10.);
    1875                 RemoveAllBoundaryPoints((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
    1876 //                double volumedifference = ConvexizeNonconvexEnvelope((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
    1877 //                double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, mol->TesselStruct, &configuration);
    1878 //                cout << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
    1879 //                cout << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
    1880                 argptr+=1;
     1899                FindNonConvexBorder((ofstream *)&cout, mol, &LCList, 5., argv[argptr+1]);
     1900//                RemoveAllBoundaryPoints((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
     1901                double volumedifference = ConvexizeNonconvexEnvelope((ofstream *)&cout, mol->TesselStruct, mol, argv[argptr]);
     1902                double clustervolume = VolumeOfConvexEnvelope((ofstream *)&cout, mol->TesselStruct, &configuration);
     1903                cout << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
     1904                cout << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
     1905                argptr+=2;
    18811906              }
    18821907              break;
  • src/moleculelist.cpp

    r8cd903 rf7f7a4  
    55 */
    66
     7#include "boundary.hpp"
    78#include "config.hpp"
    89#include "molecules.hpp"
     
    291292/** Embedding merge of a given set of molecules into one.
    292293 * Embedding merge inserts one molecule into the other.
    293  * \param *mol destination molecule
    294  * \param *srcmol source molecule
    295  * \return true - merge successful, false - merge failed (probably due to non-existant indices
    296  * \TODO find embedding center
     294 * \param *mol destination molecule (fixed one)
     295 * \param *srcmol source molecule (variable one, where atoms are taken from)
     296 * \return true - merge successful, false - merge failed (probably due to non-existant indices)
     297 * \TODO linked cell dimensions for boundary points has to be as big as inner diameter!
    297298 */
    298299bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
    299300{
    300   if (srcmol == NULL)
     301  if ((srcmol == NULL) || (mol == NULL)) {
     302    cout << Verbose(1) << "ERROR: Either fixed or variable molecule is given as NULL." << endl;
    301303    return false;
    302 
    303   // calculate center for merge
    304   srcmol->Center.CopyVector(mol->FindEmbeddingHole((ofstream *)&cout, srcmol));
    305   srcmol->Center.Zero();
    306 
    307   // perform simple merge
    308   SimpleMerge(mol, srcmol);
     304  }
     305
     306  // calculate envelope for *mol
     307  LinkedCell *LCList = new LinkedCell(mol, 8.);
     308  FindNonConvexBorder((ofstream *)&cout, mol, LCList, 4., NULL);
     309  if (mol->TesselStruct == NULL) {
     310    cout << Verbose(1) << "ERROR: Could not tesselate the fixed molecule." << endl;
     311    return false;
     312  }
     313  delete(LCList);
     314  LCList = new LinkedCell(mol->TesselStruct, 8.);  // re-create with boundary points only!
     315
     316  // prepare index list for bonds
     317  srcmol->CountAtoms((ofstream *)&cout);
     318  atom ** CopyAtoms = new atom*[srcmol->AtomCount];
     319  for(int i=0;i<srcmol->AtomCount;i++)
     320    CopyAtoms[i] = NULL;
     321
     322  // for each of the source atoms check whether we are in- or outside and add copy atom
     323  atom *Walker = srcmol->start;
     324  int nr=0;
     325  while (Walker->next != srcmol->end) {
     326    Walker = Walker->next;
     327    cout << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl;
     328    if (!mol->TesselStruct->IsInnerPoint((ofstream *)&cout, Walker->x, LCList)) {
     329      CopyAtoms[Walker->nr] = new atom(Walker);
     330      mol->AddAtom(CopyAtoms[Walker->nr]);
     331      nr++;
     332    } else {
     333      // do nothing
     334    }
     335  }
     336  cout << Verbose(1) << nr << " of " << srcmol->AtomCount << " atoms have been merged.";
     337
     338  // go through all bonds and add as well
     339  bond *Binder = srcmol->first;
     340  while(Binder->next != srcmol->last) {
     341    Binder = Binder->next;
     342    cout << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
     343    mol->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
     344  }
     345  delete(LCList);
    309346  return true;
    310347};
  • src/molecules.cpp

    r8cd903 rf7f7a4  
    662662{
    663663  int length = 0;
    664   char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
     664  const char *molname = strrchr(filename, '/');
     665  if (molname != NULL)
     666    molname += sizeof(char);  // search for filename without dirs
     667  else
     668    molname = filename; // contains no slashes
    665669  char *endname = strchr(molname, '.');
    666670  if ((endname == NULL) || (endname < molname))
     
    15621566 * \param endstep stating final configuration in molecule::Trajectories
    15631567 * \param &config configuration structure
     1568 * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
    15641569 * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
    15651570 */
    1566 bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration)
     1571bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
    15671572{
    15681573  molecule *mol = NULL;
     
    15731578  atom **PermutationMap = NULL;
    15741579  atom *Walker = NULL, *Sprinter = NULL;
    1575   MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
     1580  if (!MapByIdentity)
     1581    MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
     1582  else {
     1583    PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
     1584    Walker = start;
     1585    while (Walker->next != end) {
     1586      Walker = Walker->next;
     1587      PermutationMap[Walker->nr] = Walker;   // always pick target with the smallest distance
     1588    }
     1589  }
    15761590
    15771591  // check whether we have sufficient space in Trajectories for each atom
  • src/molecules.hpp

    r8cd903 rf7f7a4  
    182182  double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
    183183  void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
    184   bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration);
     184  bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity);
    185185       
    186186  bool CheckBounds(const Vector *x) const;
Note: See TracChangeset for help on using the changeset viewer.