source: src/mizelledata.cpp@ e2e0a5a

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/** \file mizelledata.cpp
 *
 * Implementation of a micelle creator by Daniel Dueck.
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include <iostream>
#include <fstream>
#include <config.h>
#include "atom.hpp"
#include "molecule.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/Line.hpp"
#include "World.hpp"
#include "Parser/FormatParserStorage.hpp"
#include "Parser/TremoloParser.hpp"
#include "Parser/XyzParser.hpp"
#include "Parser/PdbParser.hpp"
#include "Parser/FormatParserStorage.hpp"
#include <gsl/gsl_poly.h>
#include <gsl/gsl_eigen.h>
#include "Actions/ActionHistory.hpp"
#include "Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp"
#include "Actions/GraphAction/SubgraphDissectionAction.hpp"
#include "Shapes/BaseShapes.hpp"
#include "Shapes/ShapeOps.hpp"

#define PATH "/home/heber/tmp/"

#define AtomVector std::vector <atom *>
#define MoleculeVector std::vector <molecule *>
#define AtomList list <atom *>

int Delta2(int x1, int x2);
double Sqlength (Vector x);
int main ()
{
  setVerbosity(4);
  // need to init the history before any action is created
  ActionHistory::init();

//1.Molekuel aus pdbfile einlesen und in Variablen speichern

string path;
{
  std::ifstream file;
  path = PATH;
  path += "/tensid.data";
  file.open(path.c_str());
  FormatParserStorage::getInstance().getTremolo().load(&file);
  file.close();
}
AtomVector ever = World::getInstance().getAllAtoms();

// as all parsed atoms go into same molecule
// we don't need to create one and add them all to it
molecule *stick = ever[0]->getMolecule();

//3.Molekuel zentrieren

stick->CenterOrigin();

//4.Haupttraegheitsachse bestimmen
Vector den(0.0,0.0,1.0);

World::getInstance().clearMoleculeSelection();  // unselect all
World::getInstance().selectMolecule(stick);     // select the desired molecule for the following action
MoleculeRotateToPrincipalAxisSystem(den);
/* determine
principal axis system and make greatest eigenvector be aligned along
(0,0,1)
*/

/**/
{
  std::ofstream file;
  path = PATH;
  path += "/tensidrot.xyz";
  file.open(path.c_str());
  FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
  file.close();
}

//5.b: Molekuel um 180 Grad drehen

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

for (molecule::iterator it=stick->begin(); it !=stick->end(); ++it)
  (*it)->setPosition(RotationAxis.rotateVector((*it)->getPosition(),M_PI));



{
  std::ofstream file;
  path = PATH;
  path += "/tensid2rot.xyz";
  file.open(path.c_str());
  FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
  file.close();
}

//6.Molekuel mehrfach strukturiert mit der Haupttraegheitsachse senkrecht zu einer parametrisierten Oberflaeche anordnen

int N=200;
//6.1. Punkte auf der Oberflaeche bestimmen
//Algorithmus entnommen aus "http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere"

int ka =0;
double radius= 1.5*sqrt(pow(1.55, 2)*N);

Shape s = resize(Sphere(), radius);
std::vector<Vector> pt = s.getHomogeneousPointsOnSurface(N);

//6.2."stick" um Radius und Molekuelausdehnung in z-Richtung verschieben.

for (molecule::iterator it2=stick->begin();it2 !=stick->end();++it2) {
        Vector pos= (**it2).getPosition();
        pos[2]=pos[2]+radius; // -Min[2]
        (**it2).setPosition(pos);
}

{
  std::ofstream file;
  path = PATH;
  path += "/tensid3rot.xyz";
  file.open(path.c_str());
  FormatParserStorage::getInstance().getXyz().save(&file, World::getInstance().getAllAtoms());
  file.close();
}

//6.3.Erzeugen einer Molekuelliste, die das Molekuel "stick" "N" mal kopiert und um eine Sphaere herum verteilt 

//double MYEPSILON=1e-10;

for (ka = 0; ka<N; ka++){
  cout << "Creating " << ka+1 << " copy of tenside molecule, ";
        molecule *Tensid=stick->CopyMolecule();

        cout << "rotating ...";
  Vector ZAxis(Vector(0.0,0.0,1.0));
  Vector Axis(pt[ka]);
  const double alpha = ZAxis.Angle(Axis);
  Axis.VectorProduct(ZAxis);
  Line RotationAxis(Vector(0.0,0.0,1.0), Axis); // pt is the current Vector of point on surface

  for (molecule::iterator it2=Tensid->begin();it2 !=Tensid->end();++it2)
    (*it2)->setPosition(RotationAxis.rotateVector((*it2)->getPosition(),alpha));
  cout << "done." << endl;

  Tensid=NULL;
}

GraphSubgraphDissection();

molecule::iterator it2= stick->begin();
while(it2!=stick->end()) {
                atom *part=*it2;
                ++it2;
                stick->RemoveAtom(part);
                World::getInstance().destroyAtom (part);
  };
World::getInstance().destroyMolecule(stick);

//7. Speichern der Molekuelliste in einem data-file

{
  std::ofstream file;
  path = PATH;
  path += "/tensidoutput2.pdb";
  file.open(path.c_str());
  FormatParserStorage::getInstance().getPdb().save(&file, World::getInstance().getAllAtoms());
  file.close();
}

{
  std::ofstream file;
  path = PATH;
  path += "/tensidoutput2.data";
  file.open(path.c_str());
  FormatParserStorage::getInstance().getTremolo().save(&file, World::getInstance().getAllAtoms());
  file.close();
}

//Anhang


//gsl_matrix_free (Tensor1);
 //gsl_vector_complex_free (eval);
//gsl_matrix_complex_free (evec);
//gsl_eigen_nonsymmv_free (w);


MoleculeVector never= World::getInstance().getAllMolecules();
MoleculeVector::iterator it3 = never.begin();
while(it3!=never.end()) {
                molecule *Tensid=*it3;
                ++it3;
                World::getInstance().destroyMolecule(Tensid);
  };
ever=World::getInstance().getAllAtoms();
AtomVector::iterator it=ever.begin();
while(it!=ever.end()) {
        atom *member=*it;
        ++it;
        World::getInstance().destroyAtom(member);
        member=NULL;
}

//delete parserx;
/*for(int Zaehler=0; Zaehler<3;Zaehler++){
        delete[] Tensor[Zaehler];
        }
delete[] Tensor;
for(int Zaehler=0; Zaehler<3;Zaehler++){
        delete[] rotationmatrix[Zaehler];
}
delete[] rotationmatrix;
*/
/*for (ka = 0; ka<N; ka++){
    delete pt[ka];
}
delete[] pt;
*/

return 0;
}

int Delta2 (int x1, int x2) { 
        if (x1==x2) return 1;
        else return 0;  
}

double Sqlength (Vector x){
        double ergebnis= x[0]*x[0]+x[1]*x[1]+x[2]*x[2];
        return ergebnis;
        }

/*double[3] rotate (double[3] in[3], double theta, double phi) {
        double **rotationmatrix;
        rotationmatrix=new double *[3];
        for(int Zaehler=0; Zaehler<3;Zaehler++){
                rotationmatrix[Zaehler]=new double[3];
                }
        rotationmatrix[0][0]=cos(theta)*cos(phi);
        rotationmatrix[0][1]=-sin(phi);
        rotationmatrix[0][2]=sin(theta)*cos(phi);
        rotationmatrix[1][0]=cos(theta)*sin(phi);
        rotationmatrix[1][1]=cos(phi);
        rotationmatrix[1][2]=sin(theta)*sin(phi);
        rotationmatrix[2][0]=-sin(theta);
        rotationmatrix[2][1]=0;
        rotationmatrix[2][2]=cos(theta);
        double zahl[3];
        for(int Zaehler=0; Zaehler<3;Zaehler++){
                        Zahl[Zaehler]=rotationmatrix[Zaehler][0]*in[0]+rotationmatrix[Zaehler][1]*in[1]+rotationmatrix[Zaehler][2]*in[2];
                        }
        for(int Zaehler=0; Zaehler<3;Zaehler++){
                delete[] rotationmatrix[Zaehler];
        }
        delete[] rotationmatrix;
}*/