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(){

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