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"

namespace oopse {

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

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

DUFF::DUFF(){

    //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 = openForceFiledFile(filename);

    spMan_.parse(*ffStream, *this);

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