Ignore:
Timestamp:
Aug 6, 2008, 8:59:06 AM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
8f8621
Parents:
e8dd4d
Message:

Scaffold of new function: VerletForceIntegration() - does MD step by parsing external forces file.

PCP is too slow and unnecessarily ECut-dependent. If we incorporate the Verlet integration into molecuilder, we get a much more complete "wrapper" package.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/builder.cpp

    re8dd4d r0a08df  
    147147          third->Output(1,2,(ofstream *)&cout);
    148148          fourth->Output(1,3,(ofstream *)&cout);
    149           n.MakeNormalVector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
    150           x.CopyVector(&second->x);
     149          n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
     150          x.Copyvector(&second->x);
    151151          x.SubtractVector(&third->x);
    152           x.CopyVector(&fourth->x);
     152          x.Copyvector(&fourth->x);
    153153          x.SubtractVector(&third->x);
    154154         
    … …  
    359359      } while ((param.type = periode->FindElement(shorthand)) == NULL);
    360360      cout << Verbose(0) << "Element is " << param.type->name << endl;
    361       mol->GetAlignVector(&param);
     361      mol->GetAlignvector(&param);
    362362      for (int i=NDIM;i--;) {
    363363        x.x[i] = gsl_vector_get(param.x,i);
    … …  
    761761            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    762762            cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    763             cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
     763            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config." << endl;
    764764            cout << "\t-h/-H/-?\tGive this help screen." << endl;
    765765            cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    767767            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
    768768            cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
     769            cout << "\t-P <file> <delta_t>\tParse given forces file and append as an MD step with width delta_t to config file via Verlet." << endl;
    769770            cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
    770771            cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
    … …  
    863864              if (!mol->AddXYZFile(argv[argptr]))
    864865                cout << Verbose(2) << "File not found." << endl;
    865               else
     866              else {
    866867                cout << Verbose(2) << "File found and parsed." << endl;
    867868                config_present = present;
     869              }
    868870              break;
    869871            default:   // no match? Don't step on (this is done in next switch's default)
    … …  
    873875        if (config_present == present) {
    874876          switch(argv[argptr-1][1]) {
     877            case 'P':
     878              ExitFlag = 1;
     879              cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
     880              if (!mol->VerletForceIntegration(argv[argptr], atof(argv[argptr+1])))
     881                cout << Verbose(2) << "File not found." << endl;
     882              else
     883                cout << Verbose(2) << "File found and parsed." << endl;
     884              argptr+=2;
     885              break;
    875886            case 't':
    876887              ExitFlag = 1;
    … …  
    959970              cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
    960971              break;
     972            case 'F':
    961973            case 'f':
    962974              if (ExitFlag ==0) ExitFlag = 2;  // only set if not already by other command line switch
    … …  
    10221034                int count;
    10231035                element ** Elements;
    1024                 vector ** Vectors;
     1036                vector ** vectors;
    10251037                if (faktor < 1) {
    10261038                  cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
    … …  
    10311043                  count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
    10321044                  Elements = new element *[count];
    1033                   Vectors = new vector *[count];
     1045                  vectors = new vector *[count];
    10341046                  j = 0;
    10351047                  first = mol->start;
    … …  
    10371049                    first = first->next;
    10381050                    Elements[j] = first->type;
    1039                     Vectors[j] = &first->x;
     1051                    vectors[j] = &first->x;
    10401052                    j++;
    10411053                  }
    … …  
    10491061                    for (int k=count;k--;) { // go through every atom of the original cell
    10501062                      first = new atom(); // create a new body
    1051                       first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
     1063                      first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    10521064                      first->x.AddVector(&x);      // translate the coordinates
    10531065                      first->type = Elements[k];  // insert original element
    … …  
    10571069                  // free memory
    10581070                  delete[](Elements);
    1059                   delete[](Vectors);
     1071                  delete[](vectors);
    10601072                  // correct cell size
    10611073                  if (axis < 0) { // if sign was negative, we have to translate everything
    … …  
    11201132  clock_t start,end;
    11211133  element **Elements;
    1122   vector **Vectors;
     1134  vector **vectors;
    11231135
    11241136  // =========================== PARSE COMMAND LINE OPTIONS ====================================
    … …  
    12191231          count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
    12201232          Elements = new element *[count];
    1221           Vectors = new vector *[count];
     1233          vectors = new vector *[count];
    12221234          j = 0;
    12231235          first = mol->start;
    … …  
    12251237            first = first->next;
    12261238            Elements[j] = first->type;
    1227             Vectors[j] = &first->x;
     1239            vectors[j] = &first->x;
    12281240            j++;
    12291241          }
    … …  
    12371249            for (int k=count;k--;) { // go through every atom of the original cell
    12381250              first = new atom(); // create a new body
    1239               first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
     1251              first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    12401252              first->x.AddVector(&x);      // translate the coordinates
    12411253              first->type = Elements[k];  // insert original element
    … …  
    12471259          // free memory
    12481260          delete[](Elements);
    1249           delete[](Vectors);
     1261          delete[](vectors);
    12501262          // correct cell size
    12511263          if (axis < 0) { // if sign was negative, we have to translate everything
Note: See TracChangeset for help on using the changeset viewer.