src/UseTheForce/DUFF.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

#include "UseTheForce/DUFF.hpp"
#include "UseTheForce/DarkSide/lj_interface.h"
#include "UseTheForce/DarkSide/charge_interface.h"
#include "UseTheForce/DarkSide/dipole_interface.h"
#include "UseTheForce/DarkSide/sticky_interface.h"
#include "UseTheForce/ForceFieldFactory.hpp"
#include "io/DirectionalAtomTypesSectionParser.hpp"
#include "io/AtomTypesSectionParser.hpp"
#include "io/LennardJonesAtomTypesSectionParser.hpp"
#include "io/ElectrostaticAtomTypesSectionParser.hpp"
#include "io/EAMAtomTypesSectionParser.hpp"
#include "io/StickyAtomTypesSectionParser.hpp"
#include "io/BondTypesSectionParser.hpp"
#include "io/BendTypesSectionParser.hpp"
#include "io/TorsionTypesSectionParser.hpp"

namespace oopse {

//definition of createDUFF
ForceField* createDUFF() {
    return new DUFF();
}

//register createDUFF to ForceFieldFactory
bool registerDUFFStatus = ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);

DUFF::DUFF(){

    //set default force field filename
    setForceFieldFileName("DUFF.frc");

    //the order of adding section parsers are important
    //DirectionalAtomTypesSectionParser should be added before AtomTypesSectionParser Since
    //These two section parsers will actually create "real" AtomTypes (AtomTypesSectionParser will create
    //AtomType and DirectionalAtomTypesSectionParser will creat DirectionalAtomType which is a subclass
    //of AtomType, therefore it should come first). Other AtomTypes Section Parser will not create the 
    //"real" AtomType, they only add and set some attribute of the AtomType. Thus their order are not
    //important. AtomTypesSectionParser should be added before other atom type section parsers.
    //Make sure they are added after DirectionalAtomTypesSectionParser and AtomTypesSectionParser. 
    //The order of BondTypesSectionParser, BendTypesSectionParser and TorsionTypesSectionParser are
    //not important.
    spMan_.push_back(new DirectionalAtomTypesSectionParser());
    spMan_.push_back(new AtomTypesSectionParser());
    spMan_.push_back(new LennardJonesAtomTypesSectionParser());
    spMan_.push_back(new ElectrostaticAtomTypesSectionParser());
    spMan_.push_back(new EAMAtomTypesSectionParser());
    spMan_.push_back(new StickyAtomTypesSectionParser());
    spMan_.push_back(new BondTypesSectionParser());
    spMan_.push_back(new BendTypesSectionParser());
    spMan_.push_back(new TorsionTypesSectionParser());
    
}

void DUFF::parse(const std::string& filename) {
    ifstrstream* ffStream;
    ffStream = openForceFieldFile(filename);

    spMan_.parse(*ffStream, *this);

    ForceField::AtomTypeContainer::MapTypeIterator i;
    AtomType* at;

    for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
        at->makeFortranAtomType();
    }

    for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
        at->complete();
    }
    
}

/*
ParseState DUFF::getSection(const std::string& section) {
    ParseState result;
    
    switch(section) {
        case "AtomTypes" :
            result = DUFF::AtomTypeSection;
            break;
        case "DirectionalAtomTypes" :
            result = DUFF::DirectionalAtomTypeSection;
            break;

        case "BondTypes" :
            result = DUFF::BondTypeSection;
            break;

        case "BendTypes" :
            result = DUFF::BendTypeSection;
            break;

        case "TorsionTypes" :
            result = DUFF::TorsionTypeSection;
            break;
        default:
            result = DUFF::UnknownSection;
    }

    return result;
}

void DUFF::parse(const std::string& filename) {
    ifstrstream* ffStream;
    ffStream = openForceFieldFile(filename);
    const int bufferSize = 65535;
    std::string line;
    char buffer[bufferSize];
    int lineNo = 0;
    int atomIdent = getNAtomType() + 1;  //atom's indent begins from 1 (since only fortran's array begins from 1)
    ParseState currentSection = DUFF::UnknownSection;
    
    while(ffStream.getline(buffer, bufferSize)){
        ++lineNo;

        line = trimSpaces(buffer);
        //a line begins with "//" is comment
        if ( line.empty() || (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
            continue;
        } else {

            switch(currentSection) {
                case DUFF::AtomTypeSection :
                    parseAtomType(line, lineNo, atomIdent);
                    break;

                case DUFF::DirectionalAtomTypeSection :
                    parseDirectionalAtomType(line, lineNo);
                    break;
                    
                case DUFF::BondTypeSection :
                    parseBondType(line, lineNo);
                    break;

                case DUFF::BendTypeSection :
                    parseBendType(line, lineNo);
                    break;

                case DUFF::TorsionTypeSection :
                    parseTorsionType(line, lineNo);
                    break;
                    
                case DUFF::UnknownSection:
                    StringTokenizer tokenizer(line);
                    
                    std::string keyword = tokenizer.nextToken();
                    std::string section = tokenizer.nextToken();

                    ParseState newSection = getSection(section);
                    if (keyword != "begin" || keyword != "end") {
                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
                    } else if (keyword == "begin") {
                        if (newSection == DUFF::UnknownSection) {
                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
                        } else {
                            //enter a new section
                            currentSection = newSection;
                        }
                        
                    } else if (keyword == "end"){
                        if (currentSection == newSection) ) {
                            //leave a section
                            currentSection = DUFF::UnknownSection;
                        } else {                            
                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
                        }

                    }
                    break;
                default :
                
            }
            
        }
    }
    
    delete ffStream;
}

void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){ 
    StringTokenizer tokenizer(line);
    int nTokens = tokenizer.countTokens();    

    //in AtomTypeSection, a line at least contains 5 tokens
    //atomTypeName, is Directional, isLJ, isCharge and mass
    if (nTokens < 5)  {
                      
    } else {

        std::string atomTypeName = tokenizer.nextToken();
        bool isDirectional = tokenizer.nextTokenAsBool();
        bool isLJ = tokenizer.nextTokenAsBool();
        bool isCharge = tokenizer.nextTokenAsBool();        
        double mass = tokenizer.nextTokenAsDouble();
        double epsilon;
        double sigma;
        double charge;
        nTokens -= 5;

        //parse epsilon and sigma
        if (isLJ) {
            if (nTokens >= 2) {
                epsilon = tokenizer.nextTokenAsDouble();
                sigma = tokenizer.nextTokenAsDouble();
                nTokens -= 2;
            } else {

            }
        }

        //parse charge 
        if (isCharge) {
            if (nTokens >= 1) {
                charge = tokenizer.nextTokenAsDouble();
                nTokens -= 1;
            } else {

            }
        }
        
        AtomType* atomType;
        if (isDirectional) {
            atomType = new DirectionalAtomType();        
        } else {
            atomType = new AtomType();
        }
        
        atomType->setName(atomTypeName);
        atomType->setMass(mass);
        
        if (isLJ) {
            atomType->setLennardJones();
        }

        if (isCharge) {
            atomType->setCharge();
        }
        
        atomType->setIdent(ident); 
        
        atomType->complete();

        int setLJStatus;
        
        //notify a new LJtype atom type is created
        if (isLJ) {
            newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
        }

        int setChargeStatus;
        if (isCharge) {
            newChargeType(&ident, &charge, &setChargeStatus)
        }

        if (setLJStatus && setChargeStatus) {
            //add atom type to AtomTypeContainer
            addAtomType(atomTypeName, atomType);
            ++ident;
        } else {
            //error in notifying fortran
            delete atomType;
        }
    }    
    
}


void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
    StringTokenizer tokenizer(line);
    int nTokens = tokenizer.countTokens();  
    
    //in DirectionalAtomTypeSection, a line at least contains 6 tokens
    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
    if (nTokens < 6) {
        std::cerr << "Not enought tokens" << std::endl;
    } else {


        std::string atomTypeName = tokenizer.nextToken();
        bool isDipole = tokenizer.nextTokenAsBool();
        bool isSticky = tokenizer.nextTokenAsBool();
        double Ixx = tokenizer.nextTokenAsDouble();
        double Iyy = tokenizer.nextTokenAsDouble();
        double Izz = tokenizer.nextTokenAsDouble();
        nTokens -= 6;

        AtomType* atomType = getAtomType(atomTypeName);
        if (atomType == NULL) {

        }
        
        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
        if (dAtomType == NULL) {


        }

        if (isDipole) {
            dAtomType->setDipole();
        }

        if (isSticky) {
            dAtomType->setSticky();
        }
        
        Mat3x3d inertialMat;
        inertialMat(0, 0) = Ixx;
        inertialMat(1, 1) = Iyy;
        inertialMat(2, 2) = Izz;        
        dAtomType->setI(inertialMat);

        //read dipole moment
        double dipole;
        if (isDipole) {
            if (nTokens >= 1) {
                dipole = tokenizer.nextTokenAsDouble();
                nTokens -= 1;
            } else {

            }
        }

        //read sticky parameters
        double w0;
        double v0;
        double v0p;
        double rl;
        double ru;
        double rlp;
        double rup;        
        if (isSticky) {
            if (nTokens >= 7) {
                w0 = tokenizer.nextTokenAsDouble();
                v0 = tokenizer.nextTokenAsDouble();
                v0p = tokenizer.nextTokenAsDouble();
                rl = tokenizer.nextTokenAsDouble();
                ru = tokenizer.nextTokenAsDouble();
                rlp = tokenizer.nextTokenAsDouble();
                rup = tokenizer.nextTokenAsDouble();     
                nTokens -= 7;
            } else {

            }
        }


        //notify fotran a newDipoleType is created
        int ident = dAtomType->getIdent();
        int setDipoleStatus;
        if (isDipole) {
            newDipoleType(&ident, &dipole, &setDipoleStatus);
        }               

        //notify fotran a StickyType is created
        int setStickyStatus;
        if (isSticky) {
            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
        }


        if (!setDipoleStatus || !setStickyStatus) {

        }
        
    }
}

void DUFF::parseBondType(const std::string& line, int lineNo){ 

    StringTokenizer tokenizer(line);
    std::string at1;
    std::string at2;
    std::string bt;
    BondType* bondType = NULL;
    double b0;
    
    int nTokens = tokenizer.countTokens();

    if (nTokens < 4) {

        return;
    }
    
    at1 = tokenizer.nextToken();
    at2 = tokenizer.nextToken();
    bt = tokenizer.nextToken();
    b0 = tokenizer.nextTokenAsDouble();
    nTokens -= 4;

    //switch is a maintain nightmare
    switch(bt) {
        case "Fixed" :
            bondType = new FixedBondType();
            break;
            
        case "Harmonic" :
            if (nTokens < 1) {

            } else {

                double kb = tokenizer.nextTokenAsDouble();
                bondType = new HarmonicBondType(b0, kb);
            }

            break;

        case "Cubic" :
            if (nTokens < 4) {

            } else {

                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bondType = new CubicBondType(b0, k3, k2, k1, k0);
            }
            break;
            
        case "Quartic" :
            if (nTokens < 5) {

            } else {

                b0 = tokenizer.nextTokenAsDouble();
                double k4 = tokenizer.nextTokenAsDouble();
                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
            }
            break;

        case "Polynomial" :
            if (nTokens < 2 || nTokens % 2 != 0) {

            } else {
                int nPairs = nTokens / 2;
                int power;
                double coefficient;
                PolynomialBondType pbt = new PolynomialBondType();
                
                for (int i = 0; i < nPairs; ++i) {
                    power = tokenizer.nextTokenAsInt();
                    coefficient = tokenizer.nextTokenAsDouble();
                    pbt->setCoefficient(power, coefficient);
                }
            }
            
            break;

        default:
            
    }

    if (bondType != NULL) {
        addBondType(at1, at2, bondType);
    }
}

void DUFF::parseBendType(const std::string& line, int lineNo){ 
    StringTokenizer tokenizer(line);
    std::string at1;
    std::string at2;
    std::string at3;
    std::string bt;
    double theta0;
    BendType* bendType = NULL;

    int nTokens = tokenizer.countTokens();

    if (nTokens < 5) {

        return;
    }
    
    at1 = tokenizer.nextToken();
    at2 = tokenizer.nextToken();
    at3 = tokenizer.nextToken();
    bt = tokenizer.nextToken();
    theta0 = tokenizer.nextTokenAsDouble();
    nTokens -= 5;

    //switch is a maintain nightmare
    switch(bt) {
            
        case "Harmonic" :
            
            if (nTokens < 1) {

            } else {

                double ktheta = tokenizer.nextTokenAsDouble();
                bendType = new HarmonicBendType(theta0, ktheta);
            }
            break;
        case "GhostBend" :
            if (nTokens < 1) {

            } else {
                double ktheta = tokenizer.nextTokenAsDouble();
                bendType = new HarmonicBendType(theta0, ktheta);                
            }
            break;            

        case "UreyBradley" :
            if (nTokens < 3) {

            } else {
                double ktheta = tokenizer.nextTokenAsDouble();
                double s0 =  tokenizer.nextTokenAsDouble();
                double kub = tokenizer.nextTokenAsDouble();
                bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);                
            }
            break; 
            
        case "Cubic" :
            if (nTokens < 4) {

            } else {

                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new CubicBendType(theta0, k3, k2, k1, k0);
            }
            break;
            
        case "Quartic" :
            if (nTokens < 5) {

            } else {

                theta0 = tokenizer.nextTokenAsDouble();
                double k4 = tokenizer.nextTokenAsDouble();
                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
            }
            break;

        case "Polynomial" :
            if (nTokens < 2 || nTokens % 2 != 0) {

            } else {
                int nPairs = nTokens / 2;
                int power;
                double coefficient;
                PolynomialBendType* pbt = new PolynomialBendType();
                
                for (int i = 0; i < nPairs; ++i) {
                    power = tokenizer.nextTokenAsInt();
                    coefficient = tokenizer.nextTokenAsDouble();
                    pbt->setCoefficient(power, coefficient);
                }
            }
            
            break;

        default:
            
    }

    if (bendType != NULL) {
        addBendType(at1, at2, at3, bendType);
    }
    
}

void DUFF::parseTorsionType(const std::string& line, int lineNo){ 
    StringTokenizer tokenizer(line);
    std::string at1;
    std::string at2;
    std::string at3;
    std::string at4;
    std::string tt;
    TorsionType* torsionType = NULL;

    int nTokens = tokenizer.countTokens();

    if (nTokens < 5) {

        return;
    }
    
    at1 = tokenizer.nextToken();
    at2 = tokenizer.nextToken();
    at3 = tokenizer.nextToken();
    at4 = tokenizer.nextToken();
    tt = tokenizer.nextToken();

    nTokens -= 5;

    switch(tt) {
            
        case "Cubic" :
            if (nTokens < 4) {

            } else {

                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new CubicTorsionType(k3, k2, k1, k0);
            }
            break;
            
        case "Quartic" :
            if (nTokens < 5) {

            } else {

                theta0 = tokenizer.nextTokenAsDouble();
                double k4 = tokenizer.nextTokenAsDouble();
                double k3 = tokenizer.nextTokenAsDouble();
                double k2 = tokenizer.nextTokenAsDouble();
                double k1 = tokenizer.nextTokenAsDouble();
                double k0 = tokenizer.nextTokenAsDouble();
                
                bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
            }
            break;

        case "Polynomial" :
            if (nTokens < 2 || nTokens % 2 != 0) {

            } else {
                int nPairs = nTokens / 2;
                int power;
                double coefficient;
                PolynomialTorsionType* pbt = new PolynomialTorsionType();
                
                for (int i = 0; i < nPairs; ++i) {
                    power = tokenizer.nextTokenAsInt();
                    coefficient = tokenizer.nextTokenAsDouble();
                    pbt->setCoefficient(power, coefficient);
                }
            }
            
            break;
        case "Charmm" :
            
            if (nTokens < 3 || nTokens % 3 != 0) {

            } else {
                int nSets = nTokens / 3;
  
                CharmmTorsionType* ctt = new CharmmTorsionType();
                
                for (int i = 0; i < nSets; ++i) {
                    double kchi = tokenizer.nextTokenAsDouble();
                    int n = tokenizer.nextTokenAsInt();
                    double delta = tokenizer.nextTokenAsDouble();
    
                    ctt->setCharmmTorsionParameter(kchi, n, delta);
                }
            }
        default:
            
    }

    if (torsionType != NULL) {
        addTorsionType(at1, at2, at3, at4, torsionType);
    }
}
*/
} //end namespace oopse
Revision:	1807
Committed:	Tue Nov 30 22:43:51 2004 UTC (19 years, 9 months ago) by tim
File size:	21329 byte(s)
Log Message:	brains get built
#	User	Rev	Content
1	tim	1741	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25	gezelter	1490
26	tim	1741	#include "UseTheForce/DUFF.hpp"
27	gezelter	1634	#include "UseTheForce/DarkSide/lj_interface.h"
28	tim	1758	#include "UseTheForce/DarkSide/charge_interface.h"
29			#include "UseTheForce/DarkSide/dipole_interface.h"
30			#include "UseTheForce/DarkSide/sticky_interface.h"
31	tim	1783	#include "UseTheForce/ForceFieldFactory.hpp"
32			#include "io/DirectionalAtomTypesSectionParser.hpp"
33			#include "io/AtomTypesSectionParser.hpp"
34			#include "io/LennardJonesAtomTypesSectionParser.hpp"
35			#include "io/ElectrostaticAtomTypesSectionParser.hpp"
36			#include "io/EAMAtomTypesSectionParser.hpp"
37			#include "io/StickyAtomTypesSectionParser.hpp"
38			#include "io/BondTypesSectionParser.hpp"
39			#include "io/BendTypesSectionParser.hpp"
40			#include "io/TorsionTypesSectionParser.hpp"
41	gezelter	1490
42	tim	1741	namespace oopse {
43	gezelter	1490
44	tim	1741	//definition of createDUFF
45			ForceField* createDUFF() {
46			return new DUFF();
47			}
48	gezelter	1490
49	tim	1758	//register createDUFF to ForceFieldFactory
50	tim	1783	bool registerDUFFStatus = ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
51	gezelter	1490
52	tim	1770	DUFF::DUFF(){
53	tim	1757
54	tim	1807	//set default force field filename
55			setForceFieldFileName("DUFF.frc");
56
57	tim	1770	//the order of adding section parsers are important
58			//DirectionalAtomTypesSectionParser should be added before AtomTypesSectionParser Since
59			//These two section parsers will actually create "real" AtomTypes (AtomTypesSectionParser will create
60			//AtomType and DirectionalAtomTypesSectionParser will creat DirectionalAtomType which is a subclass
61			//of AtomType, therefore it should come first). Other AtomTypes Section Parser will not create the
62			//"real" AtomType, they only add and set some attribute of the AtomType. Thus their order are not
63			//important. AtomTypesSectionParser should be added before other atom type section parsers.
64			//Make sure they are added after DirectionalAtomTypesSectionParser and AtomTypesSectionParser.
65			//The order of BondTypesSectionParser, BendTypesSectionParser and TorsionTypesSectionParser are
66			//not important.
67			spMan_.push_back(new DirectionalAtomTypesSectionParser());
68			spMan_.push_back(new AtomTypesSectionParser());
69			spMan_.push_back(new LennardJonesAtomTypesSectionParser());
70			spMan_.push_back(new ElectrostaticAtomTypesSectionParser());
71			spMan_.push_back(new EAMAtomTypesSectionParser());
72			spMan_.push_back(new StickyAtomTypesSectionParser());
73			spMan_.push_back(new BondTypesSectionParser());
74			spMan_.push_back(new BendTypesSectionParser());
75			spMan_.push_back(new TorsionTypesSectionParser());
76
77			}
78
79			void DUFF::parse(const std::string& filename) {
80			ifstrstream* ffStream;
81	tim	1783	ffStream = openForceFieldFile(filename);
82	tim	1770
83			spMan_.parse(ffStream, this);
84
85	tim	1783	ForceField::AtomTypeContainer::MapTypeIterator i;
86	tim	1770	AtomType* at;
87
88			for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
89			at->makeFortranAtomType();
90			}
91
92			for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
93			at->complete();
94			}
95
96			}
97
98			/*
99	tim	1757	ParseState DUFF::getSection(const std::string& section) {
100			ParseState result;
101
102			switch(section) {
103			case "AtomTypes" :
104			result = DUFF::AtomTypeSection;
105			break;
106	tim	1758	case "DirectionalAtomTypes" :
107			result = DUFF::DirectionalAtomTypeSection;
108			break;
109	tim	1757
110			case "BondTypes" :
111			result = DUFF::BondTypeSection;
112			break;
113
114			case "BendTypes" :
115			result = DUFF::BendTypeSection;
116			break;
117
118			case "TorsionTypes" :
119			result = DUFF::TorsionTypeSection;
120			break;
121			default:
122			result = DUFF::UnknownSection;
123			}
124
125			return result;
126			}
127
128	tim	1741	void DUFF::parse(const std::string& filename) {
129			ifstrstream* ffStream;
130	tim	1783	ffStream = openForceFieldFile(filename);
131	tim	1741	const int bufferSize = 65535;
132	tim	1757	std::string line;
133			char buffer[bufferSize];
134	tim	1741	int lineNo = 0;
135	tim	1760	int atomIdent = getNAtomType() + 1; //atom's indent begins from 1 (since only fortran's array begins from 1)
136	tim	1757	ParseState currentSection = DUFF::UnknownSection;
137
138			while(ffStream.getline(buffer, bufferSize)){
139	tim	1741	++lineNo;
140	gezelter	1490
141	tim	1757	line = trimSpaces(buffer);
142			//a line begins with "//" is comment
143			if ( line.empty() \|\| (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
144			continue;
145			} else {
146	gezelter	1490
147	tim	1757	switch(currentSection) {
148			case DUFF::AtomTypeSection :
149	tim	1759	parseAtomType(line, lineNo, atomIdent);
150	tim	1757	break;
151	gezelter	1490
152	tim	1759	case DUFF::DirectionalAtomTypeSection :
153			parseDirectionalAtomType(line, lineNo);
154			break;
155
156	tim	1757	case DUFF::BondTypeSection :
157			parseBondType(line, lineNo);
158			break;
159	gezelter	1490
160	tim	1757	case DUFF::BendTypeSection :
161			parseBendType(line, lineNo);
162			break;
163
164			case DUFF::TorsionTypeSection :
165			parseTorsionType(line, lineNo);
166			break;
167
168			case DUFF::UnknownSection:
169			StringTokenizer tokenizer(line);
170
171			std::string keyword = tokenizer.nextToken();
172			std::string section = tokenizer.nextToken();
173
174			ParseState newSection = getSection(section);
175			if (keyword != "begin" \|\| keyword != "end") {
176			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
177			} else if (keyword == "begin") {
178			if (newSection == DUFF::UnknownSection) {
179			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
180			} else {
181			//enter a new section
182			currentSection = newSection;
183			}
184
185			} else if (keyword == "end"){
186			if (currentSection == newSection) ) {
187			//leave a section
188			currentSection = DUFF::UnknownSection;
189			} else {
190			std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
191			}
192
193			}
194			break;
195			default :
196	tim	1741
197			}
198	gezelter	1653
199	gezelter	1634	}
200	gezelter	1490	}
201
202	tim	1741	delete ffStream;
203	gezelter	1490	}
204
205	tim	1758	void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
206	tim	1757	StringTokenizer tokenizer(line);
207	tim	1758	int nTokens = tokenizer.countTokens();
208	tim	1757
209	tim	1758	//in AtomTypeSection, a line at least contains 5 tokens
210			//atomTypeName, is Directional, isLJ, isCharge and mass
211			if (nTokens < 5) {
212	tim	1759
213	tim	1758	} else {
214
215			std::string atomTypeName = tokenizer.nextToken();
216			bool isDirectional = tokenizer.nextTokenAsBool();
217			bool isLJ = tokenizer.nextTokenAsBool();
218			bool isCharge = tokenizer.nextTokenAsBool();
219			double mass = tokenizer.nextTokenAsDouble();
220			double epsilon;
221			double sigma;
222			double charge;
223			nTokens -= 5;
224
225			//parse epsilon and sigma
226			if (isLJ) {
227			if (nTokens >= 2) {
228			epsilon = tokenizer.nextTokenAsDouble();
229			sigma = tokenizer.nextTokenAsDouble();
230			nTokens -= 2;
231			} else {
232
233			}
234			}
235
236			//parse charge
237			if (isCharge) {
238			if (nTokens >= 1) {
239			charge = tokenizer.nextTokenAsDouble();
240			nTokens -= 1;
241			} else {
242
243			}
244			}
245
246			AtomType* atomType;
247			if (isDirectional) {
248			atomType = new DirectionalAtomType();
249			} else {
250			atomType = new AtomType();
251			}
252
253	tim	1757	atomType->setName(atomTypeName);
254			atomType->setMass(mass);
255
256	tim	1760	if (isLJ) {
257			atomType->setLennardJones();
258			}
259	tim	1757
260	tim	1760	if (isCharge) {
261			atomType->setCharge();
262			}
263
264	tim	1758	atomType->setIdent(ident);
265
266	tim	1757	atomType->complete();
267	tim	1758
268			int setLJStatus;
269
270	tim	1757	//notify a new LJtype atom type is created
271	tim	1758	if (isLJ) {
272			newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
273			}
274	tim	1757
275	tim	1758	int setChargeStatus;
276			if (isCharge) {
277			newChargeType(&ident, &charge, &setChargeStatus)
278			}
279
280			if (setLJStatus && setChargeStatus) {
281			//add atom type to AtomTypeContainer
282			addAtomType(atomTypeName, atomType);
283			++ident;
284			} else {
285			//error in notifying fortran
286			delete atomType;
287			}
288	tim	1757	}
289	tim	1758
290	gezelter	1490	}
291
292	tim	1758
293			void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
294			StringTokenizer tokenizer(line);
295			int nTokens = tokenizer.countTokens();
296
297	tim	1765	//in DirectionalAtomTypeSection, a line at least contains 6 tokens
298	tim	1758	//AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
299			if (nTokens < 6) {
300	tim	1765	std::cerr << "Not enought tokens" << std::endl;
301	tim	1758	} else {
302
303
304			std::string atomTypeName = tokenizer.nextToken();
305			bool isDipole = tokenizer.nextTokenAsBool();
306			bool isSticky = tokenizer.nextTokenAsBool();
307			double Ixx = tokenizer.nextTokenAsDouble();
308			double Iyy = tokenizer.nextTokenAsDouble();
309			double Izz = tokenizer.nextTokenAsDouble();
310			nTokens -= 6;
311
312			AtomType* atomType = getAtomType(atomTypeName);
313			if (atomType == NULL) {
314
315			}
316
317			DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
318			if (dAtomType == NULL) {
319
320
321			}
322
323	tim	1760	if (isDipole) {
324			dAtomType->setDipole();
325			}
326	tim	1758
327	tim	1760	if (isSticky) {
328			dAtomType->setSticky();
329			}
330
331	tim	1758	Mat3x3d inertialMat;
332			inertialMat(0, 0) = Ixx;
333	tim	1765	inertialMat(1, 1) = Iyy;
334			inertialMat(2, 2) = Izz;
335	tim	1758	dAtomType->setI(inertialMat);
336
337			//read dipole moment
338			double dipole;
339			if (isDipole) {
340			if (nTokens >= 1) {
341			dipole = tokenizer.nextTokenAsDouble();
342			nTokens -= 1;
343			} else {
344
345			}
346			}
347
348			//read sticky parameters
349			double w0;
350			double v0;
351			double v0p;
352			double rl;
353			double ru;
354			double rlp;
355			double rup;
356			if (isSticky) {
357			if (nTokens >= 7) {
358			w0 = tokenizer.nextTokenAsDouble();
359			v0 = tokenizer.nextTokenAsDouble();
360			v0p = tokenizer.nextTokenAsDouble();
361			rl = tokenizer.nextTokenAsDouble();
362			ru = tokenizer.nextTokenAsDouble();
363			rlp = tokenizer.nextTokenAsDouble();
364			rup = tokenizer.nextTokenAsDouble();
365			nTokens -= 7;
366			} else {
367
368			}
369			}
370
371
372			//notify fotran a newDipoleType is created
373			int ident = dAtomType->getIdent();
374			int setDipoleStatus;
375			if (isDipole) {
376			newDipoleType(&ident, &dipole, &setDipoleStatus);
377			}
378
379			//notify fotran a StickyType is created
380			int setStickyStatus;
381			if (isSticky) {
382			makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
383			}
384
385
386			if (!setDipoleStatus \|\| !setStickyStatus) {
387
388			}
389
390			}
391			}
392
393	tim	1757	void DUFF::parseBondType(const std::string& line, int lineNo){
394
395			StringTokenizer tokenizer(line);
396			std::string at1;
397			std::string at2;
398			std::string bt;
399			BondType* bondType = NULL;
400			double b0;
401
402			int nTokens = tokenizer.countTokens();
403
404			if (nTokens < 4) {
405
406			return;
407			}
408
409			at1 = tokenizer.nextToken();
410			at2 = tokenizer.nextToken();
411			bt = tokenizer.nextToken();
412			b0 = tokenizer.nextTokenAsDouble();
413			nTokens -= 4;
414
415			//switch is a maintain nightmare
416			switch(bt) {
417	tim	1767	case "Fixed" :
418	tim	1757	bondType = new FixedBondType();
419			break;
420
421	tim	1767	case "Harmonic" :
422	tim	1757	if (nTokens < 1) {
423
424			} else {
425
426			double kb = tokenizer.nextTokenAsDouble();
427			bondType = new HarmonicBondType(b0, kb);
428			}
429
430			break;
431
432	tim	1767	case "Cubic" :
433	tim	1757	if (nTokens < 4) {
434
435			} else {
436
437			double k3 = tokenizer.nextTokenAsDouble();
438			double k2 = tokenizer.nextTokenAsDouble();
439			double k1 = tokenizer.nextTokenAsDouble();
440			double k0 = tokenizer.nextTokenAsDouble();
441
442			bondType = new CubicBondType(b0, k3, k2, k1, k0);
443			}
444			break;
445
446	tim	1767	case "Quartic" :
447	tim	1757	if (nTokens < 5) {
448
449			} else {
450
451			b0 = tokenizer.nextTokenAsDouble();
452			double k4 = tokenizer.nextTokenAsDouble();
453			double k3 = tokenizer.nextTokenAsDouble();
454			double k2 = tokenizer.nextTokenAsDouble();
455			double k1 = tokenizer.nextTokenAsDouble();
456			double k0 = tokenizer.nextTokenAsDouble();
457
458			bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
459			}
460			break;
461
462	tim	1767	case "Polynomial" :
463	tim	1757	if (nTokens < 2 \|\| nTokens % 2 != 0) {
464
465			} else {
466			int nPairs = nTokens / 2;
467			int power;
468			double coefficient;
469			PolynomialBondType pbt = new PolynomialBondType();
470
471			for (int i = 0; i < nPairs; ++i) {
472			power = tokenizer.nextTokenAsInt();
473			coefficient = tokenizer.nextTokenAsDouble();
474	tim	1758	pbt->setCoefficient(power, coefficient);
475	tim	1757	}
476			}
477
478			break;
479
480			default:
481
482			}
483
484			if (bondType != NULL) {
485			addBondType(at1, at2, bondType);
486			}
487	gezelter	1490	}
488
489	tim	1757	void DUFF::parseBendType(const std::string& line, int lineNo){
490			StringTokenizer tokenizer(line);
491			std::string at1;
492			std::string at2;
493			std::string at3;
494			std::string bt;
495			double theta0;
496			BendType* bendType = NULL;
497
498			int nTokens = tokenizer.countTokens();
499
500			if (nTokens < 5) {
501
502			return;
503			}
504
505			at1 = tokenizer.nextToken();
506			at2 = tokenizer.nextToken();
507			at3 = tokenizer.nextToken();
508			bt = tokenizer.nextToken();
509			theta0 = tokenizer.nextTokenAsDouble();
510			nTokens -= 5;
511
512			//switch is a maintain nightmare
513			switch(bt) {
514
515	tim	1767	case "Harmonic" :
516	tim	1757
517			if (nTokens < 1) {
518
519			} else {
520
521			double ktheta = tokenizer.nextTokenAsDouble();
522			bendType = new HarmonicBendType(theta0, ktheta);
523			}
524			break;
525	tim	1767	case "GhostBend" :
526	tim	1757	if (nTokens < 1) {
527
528			} else {
529			double ktheta = tokenizer.nextTokenAsDouble();
530			bendType = new HarmonicBendType(theta0, ktheta);
531			}
532			break;
533
534	tim	1767	case "UreyBradley" :
535	tim	1757	if (nTokens < 3) {
536
537			} else {
538			double ktheta = tokenizer.nextTokenAsDouble();
539			double s0 = tokenizer.nextTokenAsDouble();
540			double kub = tokenizer.nextTokenAsDouble();
541			bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);
542			}
543			break;
544
545	tim	1767	case "Cubic" :
546	tim	1757	if (nTokens < 4) {
547
548			} else {
549
550			double k3 = tokenizer.nextTokenAsDouble();
551			double k2 = tokenizer.nextTokenAsDouble();
552			double k1 = tokenizer.nextTokenAsDouble();
553			double k0 = tokenizer.nextTokenAsDouble();
554
555			bendType = new CubicBendType(theta0, k3, k2, k1, k0);
556			}
557			break;
558
559	tim	1767	case "Quartic" :
560	tim	1757	if (nTokens < 5) {
561
562			} else {
563
564			theta0 = tokenizer.nextTokenAsDouble();
565			double k4 = tokenizer.nextTokenAsDouble();
566			double k3 = tokenizer.nextTokenAsDouble();
567			double k2 = tokenizer.nextTokenAsDouble();
568			double k1 = tokenizer.nextTokenAsDouble();
569			double k0 = tokenizer.nextTokenAsDouble();
570
571			bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
572			}
573			break;
574
575	tim	1767	case "Polynomial" :
576	tim	1757	if (nTokens < 2 \|\| nTokens % 2 != 0) {
577
578			} else {
579			int nPairs = nTokens / 2;
580			int power;
581			double coefficient;
582			PolynomialBendType* pbt = new PolynomialBendType();
583
584			for (int i = 0; i < nPairs; ++i) {
585			power = tokenizer.nextTokenAsInt();
586			coefficient = tokenizer.nextTokenAsDouble();
587	tim	1758	pbt->setCoefficient(power, coefficient);
588	tim	1757	}
589			}
590
591			break;
592
593			default:
594
595			}
596
597			if (bendType != NULL) {
598			addBendType(at1, at2, at3, bendType);
599			}
600
601	gezelter	1490	}
602
603	tim	1757	void DUFF::parseTorsionType(const std::string& line, int lineNo){
604			StringTokenizer tokenizer(line);
605			std::string at1;
606			std::string at2;
607			std::string at3;
608			std::string at4;
609			std::string tt;
610	tim	1767	TorsionType* torsionType = NULL;
611	gezelter	1490
612	tim	1757	int nTokens = tokenizer.countTokens();
613
614			if (nTokens < 5) {
615
616			return;
617			}
618
619			at1 = tokenizer.nextToken();
620			at2 = tokenizer.nextToken();
621			at3 = tokenizer.nextToken();
622			at4 = tokenizer.nextToken();
623			tt = tokenizer.nextToken();
624
625			nTokens -= 5;
626
627			switch(tt) {
628
629	tim	1767	case "Cubic" :
630	tim	1757	if (nTokens < 4) {
631
632			} else {
633
634			double k3 = tokenizer.nextTokenAsDouble();
635			double k2 = tokenizer.nextTokenAsDouble();
636			double k1 = tokenizer.nextTokenAsDouble();
637			double k0 = tokenizer.nextTokenAsDouble();
638
639			bendType = new CubicTorsionType(k3, k2, k1, k0);
640			}
641			break;
642
643	tim	1767	case "Quartic" :
644	tim	1757	if (nTokens < 5) {
645
646			} else {
647
648			theta0 = tokenizer.nextTokenAsDouble();
649			double k4 = tokenizer.nextTokenAsDouble();
650			double k3 = tokenizer.nextTokenAsDouble();
651			double k2 = tokenizer.nextTokenAsDouble();
652			double k1 = tokenizer.nextTokenAsDouble();
653			double k0 = tokenizer.nextTokenAsDouble();
654
655			bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
656			}
657			break;
658
659	tim	1767	case "Polynomial" :
660	tim	1757	if (nTokens < 2 \|\| nTokens % 2 != 0) {
661
662			} else {
663			int nPairs = nTokens / 2;
664			int power;
665			double coefficient;
666			PolynomialTorsionType* pbt = new PolynomialTorsionType();
667
668			for (int i = 0; i < nPairs; ++i) {
669			power = tokenizer.nextTokenAsInt();
670			coefficient = tokenizer.nextTokenAsDouble();
671	tim	1758	pbt->setCoefficient(power, coefficient);
672	tim	1757	}
673			}
674
675			break;
676	tim	1767	case "Charmm" :
677	tim	1757
678			if (nTokens < 3 \|\| nTokens % 3 != 0) {
679
680			} else {
681			int nSets = nTokens / 3;
682
683			CharmmTorsionType* ctt = new CharmmTorsionType();
684
685			for (int i = 0; i < nSets; ++i) {
686	tim	1758	double kchi = tokenizer.nextTokenAsDouble();
687			int n = tokenizer.nextTokenAsInt();
688			double delta = tokenizer.nextTokenAsDouble();
689	tim	1757
690	tim	1758	ctt->setCharmmTorsionParameter(kchi, n, delta);
691	tim	1757	}
692			}
693			default:
694
695			}
696
697	tim	1767	if (torsionType != NULL) {
698			addTorsionType(at1, at2, at3, at4, torsionType);
699	tim	1757	}
700	gezelter	1490	}
701	tim	1770	*/
702	tim	1741	} //end namespace oopse