ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1559 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2008-03-01 00:06:37 +00:00
parent 9346bcd722
commit 6700000007
39 changed files with 8783 additions and 8783 deletions
Download Patch File Download Diff File Expand all files Collapse all files

146
lib/poems/POEMSChain.h
View File

@ -1,73 +1,73 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: PoemsChain.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef POEMSCHAIN_H_
#define POEMSCHAIN_H_
#include "poemslist.h"
struct ChildRingData {
List<int> * childRing;
int entranceNodeId;
};
struct POEMSChain{
~POEMSChain(){
for(int i = 0; i < childChains.GetNumElements(); i++)
{
delete childChains(i);
}
}
//void printTreeStructure(int tabs);
//void getTreeAsList(List<int> * temp);
List<int> listOfNodes;
List<POEMSChain> childChains;
POEMSChain * parentChain;
List<ChildRingData> childRings;
void printTreeStructure(int tabs){
for(int i = 0; i < tabs; i++)
{
cout << "\t";
}
cout << "Chain: ";
for(int i = 0; i < listOfNodes.GetNumElements(); i++)
{
cout << *(listOfNodes(i)) << " ";
}
cout << endl;
for(int i = 0; i < childChains.GetNumElements(); i++)
{
childChains(i)->printTreeStructure(tabs + 1);
}
}
void getTreeAsList(List<int> * temp)
{
for(int i = 0; i < listOfNodes.GetNumElements(); i++)
{
int * integer = new int;
*integer = *(listOfNodes(i));
temp->Append(integer);
}
for(int i = 0; i < childChains.GetNumElements(); i++)
{
childChains(i)->getTreeAsList(temp);
}
}
};
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: PoemsChain.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef POEMSCHAIN_H_
#define POEMSCHAIN_H_
#include "poemslist.h"
struct ChildRingData {
List<int> * childRing;
int entranceNodeId;
};
struct POEMSChain{
~POEMSChain(){
for(int i = 0; i < childChains.GetNumElements(); i++)
{
delete childChains(i);
}
}
//void printTreeStructure(int tabs);
//void getTreeAsList(List<int> * temp);
List<int> listOfNodes;
List<POEMSChain> childChains;
POEMSChain * parentChain;
List<ChildRingData> childRings;
void printTreeStructure(int tabs){
for(int i = 0; i < tabs; i++)
{
cout << "\t";
}
cout << "Chain: ";
for(int i = 0; i < listOfNodes.GetNumElements(); i++)
{
cout << *(listOfNodes(i)) << " ";
}
cout << endl;
for(int i = 0; i < childChains.GetNumElements(); i++)
{
childChains(i)->printTreeStructure(tabs + 1);
}
}
void getTreeAsList(List<int> * temp)
{
for(int i = 0; i < listOfNodes.GetNumElements(); i++)
{
int * integer = new int;
*integer = *(listOfNodes(i));
temp->Append(integer);
}
for(int i = 0; i < childChains.GetNumElements(); i++)
{
childChains(i)->getTreeAsList(temp);
}
}
};
#endif

566
lib/poems/SystemProcessor.h
View File

@ -1,283 +1,283 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: SystemProcessor.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef _SYS_PROCESSOR_H_
#define _SYS_PROCESSOR_H_
#include "poemslist.h"
#include "poemstree.h"
#include "POEMSChain.h"
struct POEMSNode {
List<POEMSNode> links;
List<bool> taken;
int idNumber;
bool visited;
~POEMSNode(){
for(int i = 0; i < taken.GetNumElements(); i++)
{
delete taken(i);
}
};
};
class SystemProcessor{
private:
Tree nodes;
// List<POEMSNode> forDeletion;
List<POEMSChain> headsOfSystems;
List<List<int> > ringsInSystem;
POEMSNode * findSingleLink(TreeNode * aNode);
POEMSChain * AddNewChain(POEMSNode * currentNode);
bool setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode);
public:
SystemProcessor(void);
~SystemProcessor(void) {
headsOfSystems.DeleteValues();
for(int i = 0; i < ringsInSystem.GetNumElements(); i++)
{
for(int k = 0; k < ringsInSystem(i)->GetNumElements(); i++)
{
delete (*ringsInSystem(i))(k);
}
}
};
void processArray(int** links, int numLinks);
List<POEMSChain> * getSystemData();
int getNumberOfHeadChains();
};
SystemProcessor::SystemProcessor(void){
}
void SystemProcessor::processArray(int** links, int numLinks)
{
bool * false_var; //holds the value false; needed because a constant cannot be put into a list; the list requires a
//reference.
for(int i = 0; i < numLinks; i++) //go through all the links in the input array
{
if(!nodes.Find(links[i][0])) //if the first node in the pair is not found in the storage tree
{
POEMSNode * newNode = new POEMSNode; //make a new node
// forDeletion.Append(newNode);
newNode->idNumber = links[i][0]; //set its ID to the value
newNode->visited = false; //set it to be unvisited
nodes.Insert(links[i][0], links[i][0], (void *) newNode); //and add it to the tree storage structure
}
if(!nodes.Find(links[i][1])) //repeat process for the other half of each link
{
POEMSNode * newNode = new POEMSNode;
// forDeletion.Append(newNode);
newNode->idNumber = links[i][1];
newNode->visited = false;
nodes.Insert(links[i][1], links[i][1], (void *) newNode);
}
POEMSNode * firstNode = (POEMSNode *)nodes.Find(links[i][0]); //now that we are sure both nodes exist,
POEMSNode * secondNode = (POEMSNode *)nodes.Find(links[i][1]); //we can get both of them out of the tree
firstNode->links.Append(secondNode); //and add the link from the first to the second...
false_var = new bool;
*false_var = false; //make a new false boolean to note that the link between these two
firstNode->taken.Append(false_var); //has not already been taken, and append it to the taken list
secondNode->links.Append(firstNode); //repeat process for link from second node to first
false_var = new bool;
*false_var = false;
secondNode->taken.Append(false_var);
}
TreeNode * temp = nodes.GetRoot(); //get the root node of the node storage tree
POEMSNode * currentNode;
do
{
currentNode = findSingleLink(temp); //find the start of the next available chain
if(currentNode != NULL)
{
headsOfSystems.Append(AddNewChain(currentNode)); //and add it to the headsOfSystems list of chains
}
}
while(currentNode != NULL); //repeat this until all chains have been added
}
POEMSChain * SystemProcessor::AddNewChain(POEMSNode * currentNode){
if(currentNode == NULL) //Termination condition; if the currentNode is null, then return null
{
return NULL;
}
int * tmp;
POEMSNode * nextNode = NULL; //nextNode stores the proposed next node to add to the chain. this will be checked to make sure no backtracking is occuring before being assigned as the current node.
POEMSChain * newChain = new POEMSChain; //make a new POEMSChain object. This will be the object returned
if(currentNode->links.GetNumElements() == 0) //if we have no links from this node, then the whole chain is only one node. Add this node to the chain and return it; mark node as visited for future reference
{
currentNode->visited = true;
tmp = new int;
*tmp = currentNode->idNumber;
newChain->listOfNodes.Append(tmp);
return newChain;
}
while(currentNode->links.GetNumElements() <= 2) //we go until we get to a node that branches, or both branches have already been taken both branches can already be taken if a loop with no spurs is found in the input data
{
currentNode->visited = true;
tmp = new int;
*tmp = currentNode->idNumber;
newChain->listOfNodes.Append(tmp); //append the current node to the chain & mark as visited
//cout << "Appending node " << currentNode->idNumber << " to chain" << endl;
nextNode = currentNode->links.GetHeadElement()->value; //the next node is the first or second value stored in the links array
//of the current node. We get the first value...
if(!setLinkVisited(currentNode, nextNode)) //...and see if it points back to where we came from. If it does...
{ //either way, we set this link as visited
if(currentNode->links.GetNumElements() == 1) //if it does, then if that is the only link to this node, we're done with the chain, so append the chain to the list and return the newly created chain
{
// headsOfSystems.Append(newChain);
return newChain;
}
nextNode = currentNode->links.GetHeadElement()->next->value;//follow the other link if there is one, so we go down the chain
if(!setLinkVisited(currentNode, nextNode)) //mark link as followed, so we know not to backtrack
{
// headsOfSystems.Append(newChain);
return newChain; //This condition, where no branches have occurred but both links have already
//been taken can only occur in a loop with no spurs; add this loop to the
//system (currently added as a chain for consistency), and return.
}
}
currentNode = nextNode; //set the current node to be the next node in the chain
}
currentNode->visited = true;
tmp = new int;
*tmp = currentNode->idNumber;
newChain->listOfNodes.Append(tmp); //append the last node before branch (node shared jointly with branch chains)
//re-mark as visited, just to make sure
ListElement<POEMSNode> * tempNode = currentNode->links.GetHeadElement(); //go through all of the links, one at a time that branch
POEMSChain * tempChain = NULL; //temporary variable to hold data
while(tempNode != NULL) //when we have followed all links, stop
{
if(setLinkVisited(tempNode->value, currentNode)) //dont backtrack, or create closed loops
{
tempChain = AddNewChain(tempNode->value); //Add a new chain created out of the next node down that link
tempChain->parentChain = newChain; //set the parent to be this chain
newChain->childChains.Append(tempChain); //append the chain to this chain's list of child chains
}
tempNode = tempNode->next; //go to process the next chain
}
//headsOfSystems.Append(newChain); //append this chain to the system list
return newChain;
}
POEMSNode * SystemProcessor::findSingleLink(TreeNode * aNode)
//This function takes the root of a search tree containing POEMSNodes and returns a POEMSNode corresponding to the start of a chain in the
//system. It finds a node that has not been visited before, and only has one link; this node will be used as the head of the chain.
{
if(aNode == NULL)
{
return NULL;
}
POEMSNode * returnVal = (POEMSNode *)aNode->GetAuxData(); //get the poemsnode data out of the treenode
POEMSNode * detectLoneLoops = NULL; //is used to handle a loop that has no protruding chains
if(returnVal->visited == false)
{
detectLoneLoops = returnVal; //if we find any node that has not been visited yet, save it
}
if(returnVal->links.GetNumElements() == 1 && returnVal->visited == false) //see if it has one element and hasnt been visited already
{
return returnVal; //return the node is it meets this criteria
}
returnVal = findSingleLink(aNode->Left()); //otherwise, check the left subtree
if(returnVal == NULL) //and if we find nothing...
{
returnVal = findSingleLink(aNode->Right()); //check the right subtree
}
if(returnVal == NULL) //if we could not find any chains
{
returnVal = detectLoneLoops; //see if we found any nodes at all that havent been processed
}
return returnVal; //return what we find (will be NULL if no new chains are
//found)
}
bool SystemProcessor::setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode)
//setLinkVisited sets the links between these two nodes as visited. If they are already visited, it returns false. Otherwise, it sets
//them as visited and returns true. This function is used to see whether a certain path has been taken already in the graph structure.
//If it has been, then we need to know so we dont follow it again; this prevents infinite recursion when there is a loop, and prevents
//backtracking up a chain that has already been made. The list of booleans denoting if a link has been visited is called 'taken' and is
//part of the POEMSNode struct. The list is the same size as the list of pointers to other nodes, and stores the boolean visited/unvisited
//value for that particular link. Because each link is represented twice, (once at each node in the link), both of the boolean values need
//to be set in the event that the link has to be set as visited.
{
//cout << "Checking link between nodes " << firstNode->idNumber << " and " << secondNode->idNumber << "... ";
ListElement<POEMSNode> * tmp = firstNode->links.GetHeadElement(); //get the head element of the list of pointers for node 1
ListElement<bool> * tmp2 = firstNode->taken.GetHeadElement(); //get the head element of the list of bool isVisited flags for node 1
while(tmp->value != NULL || tmp2->value != NULL) //go through untill we reach the end of the lists
{
if(tmp->value == secondNode) //if we find the link to the other node
{
if(*(tmp2->value) == true) //if the link has already been visited
{
//cout << "visited already" << endl;
return false; //return false to indicate that the link has been visited before this attempt
}
else //otherwise, visit it
{
*tmp2->value = true;
}
break;
}
tmp = tmp->next; //go check next link
tmp2 = tmp2->next;
}
tmp = secondNode->links.GetHeadElement(); //now, if the link was unvisited, we need to go set the other node's list such that
//it also knows this link is being visited
tmp2 = secondNode->taken.GetHeadElement();
while(tmp->value != NULL || tmp2->value != NULL) //go through the list
{
if(tmp->value == firstNode) //if we find the link
{
if(*(tmp2->value) == true) //and it has already been visited, then signal an error; this shouldnt ever happen
{
cout << "Error in parsing structure! Should never reach this condition! \n" <<
"Record of visited links out of synch between two adjacent nodes.\n";
return false;
}
else
{
*tmp2->value = true; //set the appropriate value to true to indicate this link has been visited
}
break;
}
tmp = tmp->next;
tmp2 = tmp2->next;
}
//cout << "not visited" << endl;
return true; //return true to indicate that this is the first time the link has been visited
}
List<POEMSChain> * SystemProcessor::getSystemData(void) //Gets the list of POEMSChains that comprise the system. Might eventually only
//return chains linked to the reference plane, but currently returns every chain
//in the system.
{
return &headsOfSystems;
}
int SystemProcessor::getNumberOfHeadChains(void) //This function isnt implemented yet, and might be taken out entirely; this was a holdover
//from when I intended to return an array of chain pointers, rather than a list of chains
//It will probably be deleted once I finish figuring out exactly what needs to be returned
{
return 0;
}
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: SystemProcessor.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef _SYS_PROCESSOR_H_
#define _SYS_PROCESSOR_H_
#include "poemslist.h"
#include "poemstree.h"
#include "POEMSChain.h"
struct POEMSNode {
List<POEMSNode> links;
List<bool> taken;
int idNumber;
bool visited;
~POEMSNode(){
for(int i = 0; i < taken.GetNumElements(); i++)
{
delete taken(i);
}
};
};
class SystemProcessor{
private:
Tree nodes;
// List<POEMSNode> forDeletion;
List<POEMSChain> headsOfSystems;
List<List<int> > ringsInSystem;
POEMSNode * findSingleLink(TreeNode * aNode);
POEMSChain * AddNewChain(POEMSNode * currentNode);
bool setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode);
public:
SystemProcessor(void);
~SystemProcessor(void) {
headsOfSystems.DeleteValues();
for(int i = 0; i < ringsInSystem.GetNumElements(); i++)
{
for(int k = 0; k < ringsInSystem(i)->GetNumElements(); i++)
{
delete (*ringsInSystem(i))(k);
}
}
};
void processArray(int** links, int numLinks);
List<POEMSChain> * getSystemData();
int getNumberOfHeadChains();
};
SystemProcessor::SystemProcessor(void){
}
void SystemProcessor::processArray(int** links, int numLinks)
{
bool * false_var; //holds the value false; needed because a constant cannot be put into a list; the list requires a
//reference.
for(int i = 0; i < numLinks; i++) //go through all the links in the input array
{
if(!nodes.Find(links[i][0])) //if the first node in the pair is not found in the storage tree
{
POEMSNode * newNode = new POEMSNode; //make a new node
// forDeletion.Append(newNode);
newNode->idNumber = links[i][0]; //set its ID to the value
newNode->visited = false; //set it to be unvisited
nodes.Insert(links[i][0], links[i][0], (void *) newNode); //and add it to the tree storage structure
}
if(!nodes.Find(links[i][1])) //repeat process for the other half of each link
{
POEMSNode * newNode = new POEMSNode;
// forDeletion.Append(newNode);
newNode->idNumber = links[i][1];
newNode->visited = false;
nodes.Insert(links[i][1], links[i][1], (void *) newNode);
}
POEMSNode * firstNode = (POEMSNode *)nodes.Find(links[i][0]); //now that we are sure both nodes exist,
POEMSNode * secondNode = (POEMSNode *)nodes.Find(links[i][1]); //we can get both of them out of the tree
firstNode->links.Append(secondNode); //and add the link from the first to the second...
false_var = new bool;
*false_var = false; //make a new false boolean to note that the link between these two
firstNode->taken.Append(false_var); //has not already been taken, and append it to the taken list
secondNode->links.Append(firstNode); //repeat process for link from second node to first
false_var = new bool;
*false_var = false;
secondNode->taken.Append(false_var);
}
TreeNode * temp = nodes.GetRoot(); //get the root node of the node storage tree
POEMSNode * currentNode;
do
{
currentNode = findSingleLink(temp); //find the start of the next available chain
if(currentNode != NULL)
{
headsOfSystems.Append(AddNewChain(currentNode)); //and add it to the headsOfSystems list of chains
}
}
while(currentNode != NULL); //repeat this until all chains have been added
}
POEMSChain * SystemProcessor::AddNewChain(POEMSNode * currentNode){
if(currentNode == NULL) //Termination condition; if the currentNode is null, then return null
{
return NULL;
}
int * tmp;
POEMSNode * nextNode = NULL; //nextNode stores the proposed next node to add to the chain. this will be checked to make sure no backtracking is occuring before being assigned as the current node.
POEMSChain * newChain = new POEMSChain; //make a new POEMSChain object. This will be the object returned
if(currentNode->links.GetNumElements() == 0) //if we have no links from this node, then the whole chain is only one node. Add this node to the chain and return it; mark node as visited for future reference
{
currentNode->visited = true;
tmp = new int;
*tmp = currentNode->idNumber;
newChain->listOfNodes.Append(tmp);
return newChain;
}
while(currentNode->links.GetNumElements() <= 2) //we go until we get to a node that branches, or both branches have already been taken both branches can already be taken if a loop with no spurs is found in the input data
{
currentNode->visited = true;
tmp = new int;
*tmp = currentNode->idNumber;
newChain->listOfNodes.Append(tmp); //append the current node to the chain & mark as visited
//cout << "Appending node " << currentNode->idNumber << " to chain" << endl;
nextNode = currentNode->links.GetHeadElement()->value; //the next node is the first or second value stored in the links array
//of the current node. We get the first value...
if(!setLinkVisited(currentNode, nextNode)) //...and see if it points back to where we came from. If it does...
{ //either way, we set this link as visited
if(currentNode->links.GetNumElements() == 1) //if it does, then if that is the only link to this node, we're done with the chain, so append the chain to the list and return the newly created chain
{
// headsOfSystems.Append(newChain);
return newChain;
}
nextNode = currentNode->links.GetHeadElement()->next->value;//follow the other link if there is one, so we go down the chain
if(!setLinkVisited(currentNode, nextNode)) //mark link as followed, so we know not to backtrack
{
// headsOfSystems.Append(newChain);
return newChain; //This condition, where no branches have occurred but both links have already
//been taken can only occur in a loop with no spurs; add this loop to the
//system (currently added as a chain for consistency), and return.
}
}
currentNode = nextNode; //set the current node to be the next node in the chain
}
currentNode->visited = true;
tmp = new int;
*tmp = currentNode->idNumber;
newChain->listOfNodes.Append(tmp); //append the last node before branch (node shared jointly with branch chains)
//re-mark as visited, just to make sure
ListElement<POEMSNode> * tempNode = currentNode->links.GetHeadElement(); //go through all of the links, one at a time that branch
POEMSChain * tempChain = NULL; //temporary variable to hold data
while(tempNode != NULL) //when we have followed all links, stop
{
if(setLinkVisited(tempNode->value, currentNode)) //dont backtrack, or create closed loops
{
tempChain = AddNewChain(tempNode->value); //Add a new chain created out of the next node down that link
tempChain->parentChain = newChain; //set the parent to be this chain
newChain->childChains.Append(tempChain); //append the chain to this chain's list of child chains
}
tempNode = tempNode->next; //go to process the next chain
}
//headsOfSystems.Append(newChain); //append this chain to the system list
return newChain;
}
POEMSNode * SystemProcessor::findSingleLink(TreeNode * aNode)
//This function takes the root of a search tree containing POEMSNodes and returns a POEMSNode corresponding to the start of a chain in the
//system. It finds a node that has not been visited before, and only has one link; this node will be used as the head of the chain.
{
if(aNode == NULL)
{
return NULL;
}
POEMSNode * returnVal = (POEMSNode *)aNode->GetAuxData(); //get the poemsnode data out of the treenode
POEMSNode * detectLoneLoops = NULL; //is used to handle a loop that has no protruding chains
if(returnVal->visited == false)
{
detectLoneLoops = returnVal; //if we find any node that has not been visited yet, save it
}
if(returnVal->links.GetNumElements() == 1 && returnVal->visited == false) //see if it has one element and hasnt been visited already
{
return returnVal; //return the node is it meets this criteria
}
returnVal = findSingleLink(aNode->Left()); //otherwise, check the left subtree
if(returnVal == NULL) //and if we find nothing...
{
returnVal = findSingleLink(aNode->Right()); //check the right subtree
}
if(returnVal == NULL) //if we could not find any chains
{
returnVal = detectLoneLoops; //see if we found any nodes at all that havent been processed
}
return returnVal; //return what we find (will be NULL if no new chains are
//found)
}
bool SystemProcessor::setLinkVisited(POEMSNode * firstNode, POEMSNode * secondNode)
//setLinkVisited sets the links between these two nodes as visited. If they are already visited, it returns false. Otherwise, it sets
//them as visited and returns true. This function is used to see whether a certain path has been taken already in the graph structure.
//If it has been, then we need to know so we dont follow it again; this prevents infinite recursion when there is a loop, and prevents
//backtracking up a chain that has already been made. The list of booleans denoting if a link has been visited is called 'taken' and is
//part of the POEMSNode struct. The list is the same size as the list of pointers to other nodes, and stores the boolean visited/unvisited
//value for that particular link. Because each link is represented twice, (once at each node in the link), both of the boolean values need
//to be set in the event that the link has to be set as visited.
{
//cout << "Checking link between nodes " << firstNode->idNumber << " and " << secondNode->idNumber << "... ";
ListElement<POEMSNode> * tmp = firstNode->links.GetHeadElement(); //get the head element of the list of pointers for node 1
ListElement<bool> * tmp2 = firstNode->taken.GetHeadElement(); //get the head element of the list of bool isVisited flags for node 1
while(tmp->value != NULL || tmp2->value != NULL) //go through untill we reach the end of the lists
{
if(tmp->value == secondNode) //if we find the link to the other node
{
if(*(tmp2->value) == true) //if the link has already been visited
{
//cout << "visited already" << endl;
return false; //return false to indicate that the link has been visited before this attempt
}
else //otherwise, visit it
{
*tmp2->value = true;
}
break;
}
tmp = tmp->next; //go check next link
tmp2 = tmp2->next;
}
tmp = secondNode->links.GetHeadElement(); //now, if the link was unvisited, we need to go set the other node's list such that
//it also knows this link is being visited
tmp2 = secondNode->taken.GetHeadElement();
while(tmp->value != NULL || tmp2->value != NULL) //go through the list
{
if(tmp->value == firstNode) //if we find the link
{
if(*(tmp2->value) == true) //and it has already been visited, then signal an error; this shouldnt ever happen
{
cout << "Error in parsing structure! Should never reach this condition! \n" <<
"Record of visited links out of synch between two adjacent nodes.\n";
return false;
}
else
{
*tmp2->value = true; //set the appropriate value to true to indicate this link has been visited
}
break;
}
tmp = tmp->next;
tmp2 = tmp2->next;
}
//cout << "not visited" << endl;
return true; //return true to indicate that this is the first time the link has been visited
}
List<POEMSChain> * SystemProcessor::getSystemData(void) //Gets the list of POEMSChains that comprise the system. Might eventually only
//return chains linked to the reference plane, but currently returns every chain
//in the system.
{
return &headsOfSystems;
}
int SystemProcessor::getNumberOfHeadChains(void) //This function isnt implemented yet, and might be taken out entirely; this was a holdover
//from when I intended to return an array of chain pointers, rather than a list of chains
//It will probably be deleted once I finish figuring out exactly what needs to be returned
{
return 0;
}
#endif

54
lib/poems/defines.h
View File

@ -1,27 +1,27 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: defines.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef _DEFINES_H_
#define _DEFINES_H_
enum SolverType {
ONSOLVER = 0,
PARTICLESOLVER = 1
};
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: defines.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef _DEFINES_H_
#define _DEFINES_H_
enum SolverType {
ONSOLVER = 0,
PARTICLESOLVER = 1
};
#endif

2066
lib/poems/fastmatrixops.cpp
View File

File diff suppressed because it is too large Load Diff

326
lib/poems/mat4x4.cpp
View File

@ -1,163 +1,163 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: mat4x4.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include "mat4x4.h"
using namespace std;
Mat4x4::Mat4x4(){
numrows = numcols = 4;
}
Mat4x4::~Mat4x4(){
}
Mat4x4::Mat4x4(const Mat4x4& A){
numrows = numcols = 4;
elements[0][0] = A.elements[0][0];elements[0][1] = A.elements[0][1];elements[0][2] = A.elements[0][2]; elements[0][3] = A.elements[0][3];
elements[1][0] = A.elements[1][0];elements[1][1] = A.elements[1][1];elements[1][2] = A.elements[1][2]; elements[1][3] = A.elements[1][3];
elements[2][0] = A.elements[2][0];elements[2][1] = A.elements[2][1];elements[2][2] = A.elements[2][2]; elements[2][3] = A.elements[2][3];
elements[3][0] = A.elements[3][0];elements[3][1] = A.elements[3][1];elements[3][2] = A.elements[3][2]; elements[3][3] = A.elements[3][3];
}
Mat4x4::Mat4x4(const VirtualMatrix& A){
numrows = numcols = 4;
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 4) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
elements[i][j] = A.BasicGet(i,j);
}
double& Mat4x4::operator_2int(int i, int j){ // array access
return elements[i-1][j-1];
}
double Mat4x4::Get_2int(int i, int j) const{
return elements[i-1][j-1];
}
void Mat4x4::Set_2int(int i, int j, double value){
elements[i-1][j-1] = value;
}
double Mat4x4::BasicGet_2int(int i, int j) const{
return elements[i][j];
}
void Mat4x4::BasicSet_2int(int i, int j, double value){
elements[i][j] = value;
}
void Mat4x4::BasicIncrement_2int(int i, int j, double value){
elements[i][j] += value;
}
void Mat4x4::Const(double value){
elements[0][0] = elements[0][1] = elements[0][2] = elements[0][3] = value;
elements[1][0] = elements[1][1] = elements[1][2] = elements[1][3] = value;
elements[2][0] = elements[2][1] = elements[2][2] = elements[2][3] = value;
elements[3][0] = elements[3][1] = elements[3][2] = elements[3][3] = value;
}
MatrixType Mat4x4::GetType() const{
return MAT4X4;
}
istream& Mat4x4::ReadData(istream& c){ //input
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
c >> elements[i][j];
return c;
}
ostream& Mat4x4::WriteData(ostream& c) const{ //output
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
c << elements[i][j] << ' ';
return c;
}
void Mat4x4::AssignVM(const VirtualMatrix& A){
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 4) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
for(int j=0;j<numcols;j++)
elements[i][j] = A.BasicGet(i,j);
}
Mat4x4& Mat4x4::operator=(const Mat4x4& A){ // assignment operator
elements[0][0] = A.elements[0][0];elements[0][1] = A.elements[0][1];elements[0][2] = A.elements[0][2]; elements[0][3] = A.elements[0][3];
elements[1][0] = A.elements[1][0];elements[1][1] = A.elements[1][1];elements[1][2] = A.elements[1][2]; elements[1][3] = A.elements[1][3];
elements[2][0] = A.elements[2][0];elements[2][1] = A.elements[2][1];elements[2][2] = A.elements[2][2]; elements[2][3] = A.elements[2][3];
elements[3][0] = A.elements[3][0];elements[3][1] = A.elements[3][1];elements[3][2] = A.elements[3][2]; elements[3][3] = A.elements[3][3];
return *this;
}
Mat4x4& Mat4x4::operator=(const VirtualMatrix& A){ // overloaded =
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 4) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
for(int j=0;j<numcols;j++)
elements[i][j] = A.BasicGet(i,j);
return *this;
}
Mat4x4& Mat4x4::operator*=(double b){
elements[0][0] *= b; elements[0][1] *= b; elements[0][2] *= b; elements[0][3] *= b;
elements[1][0] *= b; elements[1][1] *= b; elements[1][2] *= b; elements[1][3] *= b;
elements[2][0] *= b; elements[2][1] *= b; elements[2][2] *= b; elements[2][3] *= b;
elements[3][0] *= b; elements[3][1] *= b; elements[3][2] *= b; elements[3][3] *= b;
return *this;
}
Mat4x4& Mat4x4::Identity(){
elements[0][0] = 1.0;
elements[0][1] = 0.0;
elements[0][2] = 0.0;
elements[0][3] = 0.0;
elements[1][0] = 0.0;
elements[1][1] = 1.0;
elements[1][2] = 0.0;
elements[1][3] = 0.0;
elements[2][0] = 0.0;
elements[2][1] = 0.0;
elements[2][2] = 1.0;
elements[2][3] = 0.0;
elements[3][0] = 0.0;
elements[3][1] = 0.0;
elements[3][2] = 0.0;
elements[3][3] = 1.0;
return *this;
}
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: mat4x4.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include "mat4x4.h"
using namespace std;
Mat4x4::Mat4x4(){
numrows = numcols = 4;
}
Mat4x4::~Mat4x4(){
}
Mat4x4::Mat4x4(const Mat4x4& A){
numrows = numcols = 4;
elements[0][0] = A.elements[0][0];elements[0][1] = A.elements[0][1];elements[0][2] = A.elements[0][2]; elements[0][3] = A.elements[0][3];
elements[1][0] = A.elements[1][0];elements[1][1] = A.elements[1][1];elements[1][2] = A.elements[1][2]; elements[1][3] = A.elements[1][3];
elements[2][0] = A.elements[2][0];elements[2][1] = A.elements[2][1];elements[2][2] = A.elements[2][2]; elements[2][3] = A.elements[2][3];
elements[3][0] = A.elements[3][0];elements[3][1] = A.elements[3][1];elements[3][2] = A.elements[3][2]; elements[3][3] = A.elements[3][3];
}
Mat4x4::Mat4x4(const VirtualMatrix& A){
numrows = numcols = 4;
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 4) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
elements[i][j] = A.BasicGet(i,j);
}
double& Mat4x4::operator_2int(int i, int j){ // array access
return elements[i-1][j-1];
}
double Mat4x4::Get_2int(int i, int j) const{
return elements[i-1][j-1];
}
void Mat4x4::Set_2int(int i, int j, double value){
elements[i-1][j-1] = value;
}
double Mat4x4::BasicGet_2int(int i, int j) const{
return elements[i][j];
}
void Mat4x4::BasicSet_2int(int i, int j, double value){
elements[i][j] = value;
}
void Mat4x4::BasicIncrement_2int(int i, int j, double value){
elements[i][j] += value;
}
void Mat4x4::Const(double value){
elements[0][0] = elements[0][1] = elements[0][2] = elements[0][3] = value;
elements[1][0] = elements[1][1] = elements[1][2] = elements[1][3] = value;
elements[2][0] = elements[2][1] = elements[2][2] = elements[2][3] = value;
elements[3][0] = elements[3][1] = elements[3][2] = elements[3][3] = value;
}
MatrixType Mat4x4::GetType() const{
return MAT4X4;
}
istream& Mat4x4::ReadData(istream& c){ //input
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
c >> elements[i][j];
return c;
}
ostream& Mat4x4::WriteData(ostream& c) const{ //output
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
c << elements[i][j] << ' ';
return c;
}
void Mat4x4::AssignVM(const VirtualMatrix& A){
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 4) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
for(int j=0;j<numcols;j++)
elements[i][j] = A.BasicGet(i,j);
}
Mat4x4& Mat4x4::operator=(const Mat4x4& A){ // assignment operator
elements[0][0] = A.elements[0][0];elements[0][1] = A.elements[0][1];elements[0][2] = A.elements[0][2]; elements[0][3] = A.elements[0][3];
elements[1][0] = A.elements[1][0];elements[1][1] = A.elements[1][1];elements[1][2] = A.elements[1][2]; elements[1][3] = A.elements[1][3];
elements[2][0] = A.elements[2][0];elements[2][1] = A.elements[2][1];elements[2][2] = A.elements[2][2]; elements[2][3] = A.elements[2][3];
elements[3][0] = A.elements[3][0];elements[3][1] = A.elements[3][1];elements[3][2] = A.elements[3][2]; elements[3][3] = A.elements[3][3];
return *this;
}
Mat4x4& Mat4x4::operator=(const VirtualMatrix& A){ // overloaded =
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 4) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
for(int j=0;j<numcols;j++)
elements[i][j] = A.BasicGet(i,j);
return *this;
}
Mat4x4& Mat4x4::operator*=(double b){
elements[0][0] *= b; elements[0][1] *= b; elements[0][2] *= b; elements[0][3] *= b;
elements[1][0] *= b; elements[1][1] *= b; elements[1][2] *= b; elements[1][3] *= b;
elements[2][0] *= b; elements[2][1] *= b; elements[2][2] *= b; elements[2][3] *= b;
elements[3][0] *= b; elements[3][1] *= b; elements[3][2] *= b; elements[3][3] *= b;
return *this;
}
Mat4x4& Mat4x4::Identity(){
elements[0][0] = 1.0;
elements[0][1] = 0.0;
elements[0][2] = 0.0;
elements[0][3] = 0.0;
elements[1][0] = 0.0;
elements[1][1] = 1.0;
elements[1][2] = 0.0;
elements[1][3] = 0.0;
elements[2][0] = 0.0;
elements[2][1] = 0.0;
elements[2][2] = 1.0;
elements[2][3] = 0.0;
elements[3][0] = 0.0;
elements[3][1] = 0.0;
elements[3][2] = 0.0;
elements[3][3] = 1.0;
return *this;
}

136
lib/poems/mat4x4.h
View File

@ -1,68 +1,68 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: mat4x4.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef MAT4X4_H
#define MAT4X4_H
#include "virtualmatrix.h"
#include "matrix.h"
class Vect4;
class Mat4x4 : public VirtualMatrix {
double elements[4][4];
public:
Mat4x4();
~Mat4x4();
Mat4x4(const Mat4x4& A); // copy constructor
Mat4x4(const VirtualMatrix& A); // copy constructor
double& operator_2int (int i, int j); // array access
double Get_2int(int i, int j) const;
void Set_2int(int i, int j, double value);
double BasicGet_2int(int i, int j) const;
void BasicSet_2int(int i, int j, double value);
void BasicIncrement_2int(int i, int j, double value);
void Const(double value);
MatrixType GetType() const;
std::istream& ReadData(std::istream& c); // input
std::ostream& WriteData(std::ostream& c) const; // output
void AssignVM(const VirtualMatrix& A);
Mat4x4& operator=(const Mat4x4& A); // assignment operator
Mat4x4& operator=(const VirtualMatrix& A); // overloaded =
Mat4x4& operator*=(double b);
Mat4x4& Identity();
friend Mat4x4 T(const Mat4x4& A); // a wasteful transpose
friend Mat4x4 CrossMat(Vect4& a); // a wasteful cross matrix implementation
friend void FastLU(Mat4x4& A, Mat4x4& LU, int *indx);
friend void FastLUSubs(Mat4x4& LU, Matrix& B, Matrix& C, int *indx);
friend void FastMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastMult(Mat4x4& A, Mat4x4& B, Mat4x4& C);
friend void FastMultT(Mat4x4& A, Mat4x4& B, Mat4x4& C);
friend void FastAssignT(Mat4x4& A, Mat4x4& C);
};
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: mat4x4.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef MAT4X4_H
#define MAT4X4_H
#include "virtualmatrix.h"
#include "matrix.h"
class Vect4;
class Mat4x4 : public VirtualMatrix {
double elements[4][4];
public:
Mat4x4();
~Mat4x4();
Mat4x4(const Mat4x4& A); // copy constructor
Mat4x4(const VirtualMatrix& A); // copy constructor
double& operator_2int (int i, int j); // array access
double Get_2int(int i, int j) const;
void Set_2int(int i, int j, double value);
double BasicGet_2int(int i, int j) const;
void BasicSet_2int(int i, int j, double value);
void BasicIncrement_2int(int i, int j, double value);
void Const(double value);
MatrixType GetType() const;
std::istream& ReadData(std::istream& c); // input
std::ostream& WriteData(std::ostream& c) const; // output
void AssignVM(const VirtualMatrix& A);
Mat4x4& operator=(const Mat4x4& A); // assignment operator
Mat4x4& operator=(const VirtualMatrix& A); // overloaded =
Mat4x4& operator*=(double b);
Mat4x4& Identity();
friend Mat4x4 T(const Mat4x4& A); // a wasteful transpose
friend Mat4x4 CrossMat(Vect4& a); // a wasteful cross matrix implementation
friend void FastLU(Mat4x4& A, Mat4x4& LU, int *indx);
friend void FastLUSubs(Mat4x4& LU, Matrix& B, Matrix& C, int *indx);
friend void FastMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastMult(Mat4x4& A, Mat4x4& B, Mat4x4& C);
friend void FastMultT(Mat4x4& A, Mat4x4& B, Mat4x4& C);
friend void FastAssignT(Mat4x4& A, Mat4x4& C);
};
#endif

84
lib/poems/matrices.h
View File

@ -1,42 +1,42 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: matrices.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef MATRICES_H
#define MATRICES_H
#include "matrix.h"
#include "colmatrix.h"
#include "rowmatrix.h"
#include "mat3x3.h"
#include "vect3.h"
#include "mat4x4.h"
#include "vect4.h"
#include "mat6x6.h"
#include "vect6.h"
#include "colmatmap.h"
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: matrices.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef MATRICES_H
#define MATRICES_H
#include "matrix.h"
#include "colmatrix.h"
#include "rowmatrix.h"
#include "mat3x3.h"
#include "vect3.h"
#include "mat4x4.h"
#include "vect4.h"
#include "mat6x6.h"
#include "vect6.h"
#include "colmatmap.h"
#endif

118
lib/poems/norm.cpp
View File

@ -1,59 +1,59 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: nrom.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include <cmath>
#include "norm.h"
double Magnitude(ColMatrix& A){
double G;
G = 0;
for (int i=1;i<=A.GetNumRows();i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(RowMatrix& A){
double G;
G = 0;
for (int i=1;i<=A.GetNumCols();i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(Vect3& A){
double G;
G = 0;
for (int i=1;i<=3;i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(Vect4& A){
double G;
G = 0;
for (int i=1;i<=4;i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(Vect6& A){
double G;
G = 0;
for (int i=1;i<=6;i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: nrom.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include <cmath>
#include "norm.h"
double Magnitude(ColMatrix& A){
double G;
G = 0;
for (int i=1;i<=A.GetNumRows();i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(RowMatrix& A){
double G;
G = 0;
for (int i=1;i<=A.GetNumCols();i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(Vect3& A){
double G;
G = 0;
for (int i=1;i<=3;i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(Vect4& A){
double G;
G = 0;
for (int i=1;i<=4;i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}
double Magnitude(Vect6& A){
double G;
G = 0;
for (int i=1;i<=6;i++) G += A.Get(i)*A.Get(i);
G = sqrt(G);
return G;
}

60
lib/poems/norm.h
View File

@ -1,30 +1,30 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: norm.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef NORM_H
#define NORM_H
#include "matrices.h"
double Magnitude(ColMatrix& A);
double Magnitude(RowMatrix& A);
double Magnitude(Vect3& A);
double Magnitude(Vect4& A);
double Magnitude(Vect6& A);
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: norm.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef NORM_H
#define NORM_H
#include "matrices.h"
double Magnitude(ColMatrix& A);
double Magnitude(RowMatrix& A);
double Magnitude(Vect3& A);
double Magnitude(Vect4& A);
double Magnitude(Vect6& A);
#endif

324
lib/poems/poemsnodelib.h
View File

@ -1,162 +1,162 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: poemsnodelib.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef NODELIB_H
#define NODELIB_H
#include <iostream>
using namespace std;
TreeNode *GetTreeNode(int item,TreeNode *lptr = NULL,TreeNode *rptr =NULL);
void FreeTreeNode(TreeNode *p);
void Postorder(TreeNode *t, void visit(TreeNode* &t));
void Preorder (TreeNode *t, void visit(TreeNode* &t));
void CountLeaf (TreeNode *t, int& count);
int Depth (TreeNode *t);
void IndentBlanks(int num);
void PrintTree (TreeNode *t, int level);
// ---------------- Global functions-----------------//
// postorder recursive scan of the nodes in a tree
void Postorder (TreeNode *t, void visit(TreeNode* &t))
{
// the recursive scan terminates on a empty subtree
if (t != NULL)
{
Postorder(t->Left(), visit); // descend left
Postorder(t->Right(), visit); // descend right
visit(t); // visit the node
}
}
// preorder recursive scan of the nodes in a tree
void Preorder (TreeNode *t, void visit(TreeNode* &t))
{
// the recursive scan terminates on a empty subtree
if (t != NULL)
{
visit(t); // visit the node
Preorder(t->Left(), visit); // descend left
Preorder(t->Right(), visit); // descend right
}
}
//create TreeNode object with pointer fields lptr and rptr
// The pointers have default value NULL
TreeNode *GetTreeNode(int item,TreeNode *lptr,TreeNode *rptr)
{
TreeNode *p;
// call new to allocate the new node
// pass parameters lptr and rptr to the function
p = new TreeNode(item, lptr, rptr);
// if insufficient memory, terminatewith an error message
if (p == NULL)
{
cerr << "Memory allocation failure!\n";
exit(1);
}
// return the pointer to the system generated memory
return p;
}
// deallocate dynamic memory associated with the node
void FreeTreeNode(TreeNode *p)
{
delete p;
}
// the function uses the postorder scan. a visit
// tests whether the node is a leaf node
void CountLeaf (TreeNode *t, int& count)
{
//use postorder descent
if(t !=NULL)
{
CountLeaf(t->Left(), count); // descend left
CountLeaf(t->Right(), count); // descend right
// check if node t is a leaf node (no descendants)
// if so, increment the variable count
if (t->Left() == NULL && t->Right() == NULL)
count++;
}
}
// the function uses the postorder scan. it computes the
// depth of the left and right subtrees of a node and
// returns the depth as 1 + max(depthLeft,depthRight).
// the depth of an empty tree is -1
int Depth (TreeNode *t)
{
int depthLeft, depthRight, depthval;
if (t == NULL)
depthval = -1;
else
{
depthLeft = Depth(t->Left());
depthRight = Depth(t->Right());
depthval = 1+(depthLeft > depthRight?depthLeft:depthRight);
}
return depthval;
}
void IndentBlanks(int num)
{
// const int indentblock = 6;
for(int i = 0; i < num; i++)
cout << " ";
}
void PrintTree (TreeNode *t, int level)
{
//print tree with root t, as long as t!=NULL
if (t != NULL)
{
int indentUnit = 5;
// print right branch of tree t
PrintTree(t->Right(),level + 1);
// indent to current level; output node data
IndentBlanks(indentUnit*level);
cout << t->GetData() << endl;
// print left branch of tree t
PrintTree(t->Left(),level + 1);
}
}
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: poemsnodelib.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef NODELIB_H
#define NODELIB_H
#include <iostream>
using namespace std;
TreeNode *GetTreeNode(int item,TreeNode *lptr = NULL,TreeNode *rptr =NULL);
void FreeTreeNode(TreeNode *p);
void Postorder(TreeNode *t, void visit(TreeNode* &t));
void Preorder (TreeNode *t, void visit(TreeNode* &t));
void CountLeaf (TreeNode *t, int& count);
int Depth (TreeNode *t);
void IndentBlanks(int num);
void PrintTree (TreeNode *t, int level);
// ---------------- Global functions-----------------//
// postorder recursive scan of the nodes in a tree
void Postorder (TreeNode *t, void visit(TreeNode* &t))
{
// the recursive scan terminates on a empty subtree
if (t != NULL)
{
Postorder(t->Left(), visit); // descend left
Postorder(t->Right(), visit); // descend right
visit(t); // visit the node
}
}
// preorder recursive scan of the nodes in a tree
void Preorder (TreeNode *t, void visit(TreeNode* &t))
{
// the recursive scan terminates on a empty subtree
if (t != NULL)
{
visit(t); // visit the node
Preorder(t->Left(), visit); // descend left
Preorder(t->Right(), visit); // descend right
}
}
//create TreeNode object with pointer fields lptr and rptr
// The pointers have default value NULL
TreeNode *GetTreeNode(int item,TreeNode *lptr,TreeNode *rptr)
{
TreeNode *p;
// call new to allocate the new node
// pass parameters lptr and rptr to the function
p = new TreeNode(item, lptr, rptr);
// if insufficient memory, terminatewith an error message
if (p == NULL)
{
cerr << "Memory allocation failure!\n";
exit(1);
}
// return the pointer to the system generated memory
return p;
}
// deallocate dynamic memory associated with the node
void FreeTreeNode(TreeNode *p)
{
delete p;
}
// the function uses the postorder scan. a visit
// tests whether the node is a leaf node
void CountLeaf (TreeNode *t, int& count)
{
//use postorder descent
if(t !=NULL)
{
CountLeaf(t->Left(), count); // descend left
CountLeaf(t->Right(), count); // descend right
// check if node t is a leaf node (no descendants)
// if so, increment the variable count
if (t->Left() == NULL && t->Right() == NULL)
count++;
}
}
// the function uses the postorder scan. it computes the
// depth of the left and right subtrees of a node and
// returns the depth as 1 + max(depthLeft,depthRight).
// the depth of an empty tree is -1
int Depth (TreeNode *t)
{
int depthLeft, depthRight, depthval;
if (t == NULL)
depthval = -1;
else
{
depthLeft = Depth(t->Left());
depthRight = Depth(t->Right());
depthval = 1+(depthLeft > depthRight?depthLeft:depthRight);
}
return depthval;
}
void IndentBlanks(int num)
{
// const int indentblock = 6;
for(int i = 0; i < num; i++)
cout << " ";
}
void PrintTree (TreeNode *t, int level)
{
//print tree with root t, as long as t!=NULL
if (t != NULL)
{
int indentUnit = 5;
// print right branch of tree t
PrintTree(t->Right(),level + 1);
// indent to current level; output node data
IndentBlanks(indentUnit*level);
cout << t->GetData() << endl;
// print left branch of tree t
PrintTree(t->Left(),level + 1);
}
}
#endif

270
lib/poems/vect4.cpp
View File

@ -1,135 +1,135 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: vect4.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include "vect4.h"
using namespace std;
Vect4::Vect4(){
numrows = 4; numcols = 1;
}
Vect4::~Vect4(){
}
Vect4::Vect4(const Vect4& A){ // copy constructor
numrows = 4; numcols = 1;
elements[0] = A.elements[0];
elements[1] = A.elements[1];
elements[2] = A.elements[2];
elements[3] = A.elements[3];
}
Vect4::Vect4(const VirtualMatrix& A){ // copy constructor
numrows = 4; numcols = 1;
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 1) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<4;i++)
elements[i] = A.BasicGet(i,0);
}
double& Vect4::operator_1int (int i){ // array access
return elements[i-1];
}
double Vect4::Get_1int(int i) const{
return elements[i-1];
}
void Vect4::Set_1int(int i, double value){
elements[i-1] = value;
}
double Vect4::BasicGet_1int(int i) const{
return elements[i];
}
void Vect4::BasicSet_1int(int i, double value){
elements[i] = value;
}
void Vect4::BasicIncrement_1int(int i, double value){
elements[i] += value;
}
void Vect4::Const(double value){
elements[0] = value;
elements[1] = value;
elements[2] = value;
elements[3] = value;
}
MatrixType Vect4::GetType() const{
return VECT4;
}
istream& Vect4::ReadData(istream& c){ //input
for(int i=0;i<4;i++)
c >> elements[i];
return c;
}
ostream& Vect4::WriteData(ostream& c) const{ //output
for(int i=0;i<4;i++)
c << elements[i] << ' ';
return c;
}
void Vect4::AssignVM(const VirtualMatrix& A){
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 1) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
elements[i] = A.BasicGet(i,0);
}
Vect4& Vect4::operator=(const Vect4& A){ // assignment operator
elements[0] = A.elements[0];
elements[1] = A.elements[1];
elements[2] = A.elements[2];
elements[3] = A.elements[3];
return *this;
}
Vect4& Vect4::operator=(const VirtualMatrix& A){ // overloaded =
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 1) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
elements[i] = A.BasicGet(i,0);
return *this;
}
Vect4& Vect4::operator*=(double b){
elements[0] *= b;
elements[1] *= b;
elements[2] *= b;
elements[3] *= b;
return *this;
}
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: vect4.cpp *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#include "vect4.h"
using namespace std;
Vect4::Vect4(){
numrows = 4; numcols = 1;
}
Vect4::~Vect4(){
}
Vect4::Vect4(const Vect4& A){ // copy constructor
numrows = 4; numcols = 1;
elements[0] = A.elements[0];
elements[1] = A.elements[1];
elements[2] = A.elements[2];
elements[3] = A.elements[3];
}
Vect4::Vect4(const VirtualMatrix& A){ // copy constructor
numrows = 4; numcols = 1;
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 1) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<4;i++)
elements[i] = A.BasicGet(i,0);
}
double& Vect4::operator_1int (int i){ // array access
return elements[i-1];
}
double Vect4::Get_1int(int i) const{
return elements[i-1];
}
void Vect4::Set_1int(int i, double value){
elements[i-1] = value;
}
double Vect4::BasicGet_1int(int i) const{
return elements[i];
}
void Vect4::BasicSet_1int(int i, double value){
elements[i] = value;
}
void Vect4::BasicIncrement_1int(int i, double value){
elements[i] += value;
}
void Vect4::Const(double value){
elements[0] = value;
elements[1] = value;
elements[2] = value;
elements[3] = value;
}
MatrixType Vect4::GetType() const{
return VECT4;
}
istream& Vect4::ReadData(istream& c){ //input
for(int i=0;i<4;i++)
c >> elements[i];
return c;
}
ostream& Vect4::WriteData(ostream& c) const{ //output
for(int i=0;i<4;i++)
c << elements[i] << ' ';
return c;
}
void Vect4::AssignVM(const VirtualMatrix& A){
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 1) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
elements[i] = A.BasicGet(i,0);
}
Vect4& Vect4::operator=(const Vect4& A){ // assignment operator
elements[0] = A.elements[0];
elements[1] = A.elements[1];
elements[2] = A.elements[2];
elements[3] = A.elements[3];
return *this;
}
Vect4& Vect4::operator=(const VirtualMatrix& A){ // overloaded =
// error check
if( (A.GetNumRows() != 4) || (A.GetNumCols() != 1) ){
cerr << "illegal matrix size" << endl;
exit(0);
}
for(int i=0;i<numrows;i++)
elements[i] = A.BasicGet(i,0);
return *this;
}
Vect4& Vect4::operator*=(double b){
elements[0] *= b;
elements[1] *= b;
elements[2] *= b;
elements[3] *= b;
return *this;
}

142
lib/poems/vect4.h
View File

@ -1,71 +1,71 @@
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: vect4.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef VECT4_H
#define VECT4_H
#include "virtualcolmatrix.h"
class Matrix;
class Mat4x4;
class Vect4 : public VirtualColMatrix {
double elements[4];
public:
Vect4();
~Vect4();
Vect4(const Vect4& A); // copy constructor
Vect4(const VirtualMatrix& A); // copy constructor
double& operator_1int(int i); // array access
double Get_1int(int i) const;
void Set_1int(int i, double value);
double BasicGet_1int(int i) const;
void BasicSet_1int(int i, double value);
void BasicIncrement_1int(int i, double value);
void Const(double value);
MatrixType GetType() const;
std::ostream& WriteData(std::ostream& c) const;
std::istream& ReadData(std::istream& c);
void AssignVM(const VirtualMatrix& A);
Vect4& operator=(const Vect4& A); // assignment operator
Vect4& operator=(const VirtualMatrix& A); // overloaded =
Vect4& operator*=(double b);
friend Matrix T(const Vect4& A); // a wasteful transpose
friend Mat4x4 CrossMat(Vect4& a); // a wasteful cross matrix implementation
// fast matrix functions
friend void FastAssign(Vect4& a, Vect4& c);
friend void FastSimpleRotation(Vect4& v, double q, Mat4x4& d);
friend void FastCross(Vect4& a, Vect4& b, Vect4& c); // cross product axb = c
friend void FastTripleSum(Vect4& a, Vect4& b, Vect4& c, Vect4& d);
friend void FastTripleSumPPM(Vect4& a, Vect4& b, Vect4& c, Vect4& d);
friend void FastMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastMult(double a, Vect4& B, Vect4& C);
friend void FastAdd(Vect4& A, Vect4& B, Vect4& C);
friend void FastSubt(Vect4& A, Vect4& B, Vect4& C);
};
#endif
/*
*_________________________________________________________________________*
* POEMS: PARALLELIZABLE OPEN SOURCE EFFICIENT MULTIBODY SOFTWARE *
* DESCRIPTION: SEE READ-ME *
* FILE NAME: vect4.h *
* AUTHORS: See Author List *
* GRANTS: See Grants List *
* COPYRIGHT: (C) 2005 by Authors as listed in Author's List *
* LICENSE: Please see License Agreement *
* DOWNLOAD: Free at www.rpi.edu/~anderk5 *
* ADMINISTRATOR: Prof. Kurt Anderson *
* Computational Dynamics Lab *
* Rensselaer Polytechnic Institute *
* 110 8th St. Troy NY 12180 *
* CONTACT: anderk5@rpi.edu *
*_________________________________________________________________________*/
#ifndef VECT4_H
#define VECT4_H
#include "virtualcolmatrix.h"
class Matrix;
class Mat4x4;
class Vect4 : public VirtualColMatrix {
double elements[4];
public:
Vect4();
~Vect4();
Vect4(const Vect4& A); // copy constructor
Vect4(const VirtualMatrix& A); // copy constructor
double& operator_1int(int i); // array access
double Get_1int(int i) const;
void Set_1int(int i, double value);
double BasicGet_1int(int i) const;
void BasicSet_1int(int i, double value);
void BasicIncrement_1int(int i, double value);
void Const(double value);
MatrixType GetType() const;
std::ostream& WriteData(std::ostream& c) const;
std::istream& ReadData(std::istream& c);
void AssignVM(const VirtualMatrix& A);
Vect4& operator=(const Vect4& A); // assignment operator
Vect4& operator=(const VirtualMatrix& A); // overloaded =
Vect4& operator*=(double b);
friend Matrix T(const Vect4& A); // a wasteful transpose
friend Mat4x4 CrossMat(Vect4& a); // a wasteful cross matrix implementation
// fast matrix functions
friend void FastAssign(Vect4& a, Vect4& c);
friend void FastSimpleRotation(Vect4& v, double q, Mat4x4& d);
friend void FastCross(Vect4& a, Vect4& b, Vect4& c); // cross product axb = c
friend void FastTripleSum(Vect4& a, Vect4& b, Vect4& c, Vect4& d);
friend void FastTripleSumPPM(Vect4& a, Vect4& b, Vect4& c, Vect4& d);
friend void FastMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastNegTMult(Mat4x4& A, Vect4& B, Vect4& C);
friend void FastMult(double a, Vect4& B, Vect4& C);
friend void FastAdd(Vect4& A, Vect4& B, Vect4& C);
friend void FastSubt(Vect4& A, Vect4& B, Vect4& C);
};
#endif

44
potentials/GaN.sw
View File

@ -1,22 +1,22 @@
# GaN potential: A. Bere, and A. Serra, Phil. Mag. 86, 2159(2006)
# note that the parameters for this literature potential are pairwise
# so that there are some flexibility in the way the
# parameters can be entered. As one way, we assume that
# lambda_ijk is equal to lambda_ik and eps_ijk is
# equal to sqrt(lambda_ij*eps_ij*lambda_ik*eps_ik)/lambda_ik,
# and all other parameters in the ijk line are for ik.
# These entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
#
# eps sigma a lambda gamma cos(theta) A B p q tol
#
Ga Ga Ga 1.2000000 2.100 1.6 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N N N 1.2000000 1.300 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
Ga Ga N 1.6136914 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
Ga N N 2.1700000 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N Ga Ga 2.1700000 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N Ga N 1.6136914 1.300 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N N Ga 1.6136914 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
Ga N Ga 1.6136914 2.100 1.6 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
# GaN potential: A. Bere, and A. Serra, Phil. Mag. 86, 2159(2006)
# note that the parameters for this literature potential are pairwise
# so that there are some flexibility in the way the
# parameters can be entered. As one way, we assume that
# lambda_ijk is equal to lambda_ik and eps_ijk is
# equal to sqrt(lambda_ij*eps_ij*lambda_ik*eps_ik)/lambda_ik,
# and all other parameters in the ijk line are for ik.
# These entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
#
# eps sigma a lambda gamma cos(theta) A B p q tol
#
Ga Ga Ga 1.2000000 2.100 1.6 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N N N 1.2000000 1.300 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
Ga Ga N 1.6136914 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
Ga N N 2.1700000 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N Ga Ga 2.1700000 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N Ga N 1.6136914 1.300 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
N N Ga 1.6136914 1.695 1.8 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0
Ga N Ga 1.6136914 2.100 1.6 32.5 1.2 -0.333333333333 7.917 0.72 4.0 0.0 0.0

34
potentials/Si.sw
View File

@ -1,17 +1,17 @@
# Stillinger-Weber parameters for various elements and mixtures
# multiple entries can be added to this file, LAMMPS reads the ones it needs
# these entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
# format of a single entry (one or more lines):
# element 1, element 2, element 3,
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
# Here are the original parameters in metal units, for Silicon from:
#
# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985)
#
Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333
7.049556277 0.6022245584 4.0 0.0 0.0
# Stillinger-Weber parameters for various elements and mixtures
# multiple entries can be added to this file, LAMMPS reads the ones it needs
# these entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
# format of a single entry (one or more lines):
# element 1, element 2, element 3,
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
# Here are the original parameters in metal units, for Silicon from:
#
# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985)
#
Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333
7.049556277 0.6022245584 4.0 0.0 0.0

100
src/MAKE/Windows/erfc.cpp
View File

@ -1,51 +1,51 @@
//This code was written by Philip Nicoletti
//http://www.codeguru.com/forum/archive/index.php/t-129990.html
//
//Modified by Jin Ma, Oklahoma State University for LAMMPS
//erfc() is defined in GNU libraries. This code is a simplified
//version for implementation with Visual C++.
//
//Warning: these functions are not fully tested.
//
#include "erfc.h"
#include "math.h"
double erf(double x)
{
//
// Computation of the error function erf(x).
//
return (1-erfc(x));
}
//
//
double erfc(double x)
{
//
// Computation of the complementary error function erfc(x).
//
// The algorithm is based on a Chebyshev fit as denoted in
// Numerical Recipes 2nd ed. on p. 214 (W.H.Press et al.).
//
// The fractional error is always less than 1.2e-7.
//
//
// The parameters of the Chebyshev fit
//
const double a1 = -1.26551223, a2 = 1.00002368,
a3 = 0.37409196, a4 = 0.09678418,
a5 = -0.18628806, a6 = 0.27886807,
a7 = -1.13520398, a8 = 1.48851587,
a9 = -0.82215223, a10 = 0.17087277;
//
double v = 1; // The return value
double z = fabs(x);
//
if (z == 0) return v; // erfc(0)=1
double t = 1/(1+0.5*z);
v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+
t*(a7+t*(a8+t*(a9+t*a10)))))))));
if (x < 0) v = 2-v; // erfc(-x)=2-erfc(x)
return v;
//This code was written by Philip Nicoletti
//http://www.codeguru.com/forum/archive/index.php/t-129990.html
//
//Modified by Jin Ma, Oklahoma State University for LAMMPS
//erfc() is defined in GNU libraries. This code is a simplified
//version for implementation with Visual C++.
//
//Warning: these functions are not fully tested.
//
#include "erfc.h"
#include "math.h"
double erf(double x)
{
//
// Computation of the error function erf(x).
//
return (1-erfc(x));
}
//
//
double erfc(double x)
{
//
// Computation of the complementary error function erfc(x).
//
// The algorithm is based on a Chebyshev fit as denoted in
// Numerical Recipes 2nd ed. on p. 214 (W.H.Press et al.).
//
// The fractional error is always less than 1.2e-7.
//
//
// The parameters of the Chebyshev fit
//
const double a1 = -1.26551223, a2 = 1.00002368,
a3 = 0.37409196, a4 = 0.09678418,
a5 = -0.18628806, a6 = 0.27886807,
a7 = -1.13520398, a8 = 1.48851587,
a9 = -0.82215223, a10 = 0.17087277;
//
double v = 1; // The return value
double z = fabs(x);
//
if (z == 0) return v; // erfc(0)=1
double t = 1/(1+0.5*z);
v = t*exp((-z*z) +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+
t*(a7+t*(a8+t*(a9+t*a10)))))))));
if (x < 0) v = 2-v; // erfc(-x)=2-erfc(x)
return v;
}

8
src/MAKE/Windows/erfc.h
View File

@ -1,5 +1,5 @@
//
double erf(double x);
//
double erf(double x);
double erfc(double x);

252
src/MAKE/Windows/notes.2
View File

@ -1,126 +1,126 @@
This is instruction for the modification of LAMMPS for MS Windows
LAMMPS version: Feb 2007
compiled without MPI and FFT in Viusal C++ 6.0
(All packages except for XTC, MEAM appear to work.)
-------------------
1. Create an empty workspace (Win32 console), add all .h and .cpp
files into the project.
2. At about 80 places in the code, variables are redefined. Most of
these variables are loop counters, which can be easily fixed.
Code looks like this:
for (int i=0; i<5; i++) {
something;
}
for (int i=0; i<5; i++) {
something else;
}
This is ok with g++ compiler. But VC thinks the i is redefined in the
second loop. So the variable scope is different. This happens many times
in the code. It can be fixed easily based on the compiling error.
3. At the beginning of fft3d.h, added:
#ifndef FFT_NONE
#define FFT_NONE
#endif
4. In input.cpp, changed the two header files
//#include "unistd.h"
#include "direct.h"
4A. (added by Tim Lau, MIT, ttl@mit.edu)
In variable.cpp, change the header files
//#include "unistd.h"
#include "sleep.h"
Add in the included sleep.h and sleep.cpp files.
4B. (added by Tim Lau, MIT, ttl@mit.edu)
In shell.cpp, change the header file:
//#include "unistd.h"
#include "direct.h"
Change the line in shell.cpp:
mkdir(arg[i], S_IRWXU | S_IRGRP | S_IXGRP);
to:
mkdir(arg[i]);
since Windows obviously does not use UNIX file permissions.
It's also possible that the line has to be changed to:
_mkdir(arg[i]);
depending on the version of the Visual C++ compiler used.
5. Added mpi.h and mpi.cpp (in STUBS folder) to the workspace
In mpi.cpp, commented the time.h header file
//#include <sys/time.h>
commented the original code in MPI_Wtime(), just make it return 0;
6. In system.cpp, two changes due to difference in the input argument
list
Line 83: int iarg = 2;
Line 172: inflag=1; //add this line
The number of input arguments (nargs) is different in g++ and VC when
you give arguments to run a program. This might be related to MPI as
well. The difference is one. Once the above changes are made, the
program is taking the correct argument.
However, it has been observed in the latest versions of sytem.cpp that
no modification needs be made to the file as distributed from the
LAMMPS website to work. The user however, instead of starting LAMMPS
by the command:
lammps in.file
as he would if he implemented the changes detailed here, would launch
in the Unix style:
lammps < in.file
7. The new version LAMMPS calls the error function:
double erfc(double)
This function is in the GNU C library. However, it's not found for
VC++.
Three options:
a. One can try to find erfc() from other libraries.
b. The erfc() is called for pair_modify table option. One can set
the table option to be 0 to avoid calling this function.
c. Write your own functions.
In this code, two files erfc.h, erfc.cpp are created and added to the project.
Files that call erfc() all add
#include "erfc.h" at the beginning.
Note: the functions are not fully tested, use with caution.
8. MSVC does not have a inttypes.h file. The simplest way
to deal with this problem is to download inttypes.h from the
following site:
http://www.koders.com/c/fidDE7D6EFFD475FAB1B7F6A2BBA791401CFA88FFA3.aspx
and add this file into the workspace.
9. MSVC does not have dirent.h. The problem is solved by downloading
a version of it for Windows from the following website:
http://www.softagalleria.net/dirent/index.en.html
10. Every time an error pops up for a line like this:
char *words[params_per_line];
replace to the following type:
char **words = new char*[params_per_line];
The dynamic memory allocation in MSVC requires a line like this.
11. Build the project. Specify appropriate input file to run the code.
The Windows result might be different from Unix results. Be Cautious.
This is instruction for the modification of LAMMPS for MS Windows
LAMMPS version: Feb 2007
compiled without MPI and FFT in Viusal C++ 6.0
(All packages except for XTC, MEAM appear to work.)
-------------------
1. Create an empty workspace (Win32 console), add all .h and .cpp
files into the project.
2. At about 80 places in the code, variables are redefined. Most of
these variables are loop counters, which can be easily fixed.
Code looks like this:
for (int i=0; i<5; i++) {
something;
}
for (int i=0; i<5; i++) {
something else;
}
This is ok with g++ compiler. But VC thinks the i is redefined in the
second loop. So the variable scope is different. This happens many times
in the code. It can be fixed easily based on the compiling error.
3. At the beginning of fft3d.h, added:
#ifndef FFT_NONE
#define FFT_NONE
#endif
4. In input.cpp, changed the two header files
//#include "unistd.h"
#include "direct.h"
4A. (added by Tim Lau, MIT, ttl@mit.edu)
In variable.cpp, change the header files
//#include "unistd.h"
#include "sleep.h"
Add in the included sleep.h and sleep.cpp files.
4B. (added by Tim Lau, MIT, ttl@mit.edu)
In shell.cpp, change the header file:
//#include "unistd.h"
#include "direct.h"
Change the line in shell.cpp:
mkdir(arg[i], S_IRWXU | S_IRGRP | S_IXGRP);
to:
mkdir(arg[i]);
since Windows obviously does not use UNIX file permissions.
It's also possible that the line has to be changed to:
_mkdir(arg[i]);
depending on the version of the Visual C++ compiler used.
5. Added mpi.h and mpi.cpp (in STUBS folder) to the workspace
In mpi.cpp, commented the time.h header file
//#include <sys/time.h>
commented the original code in MPI_Wtime(), just make it return 0;
6. In system.cpp, two changes due to difference in the input argument
list
Line 83: int iarg = 2;
Line 172: inflag=1; //add this line
The number of input arguments (nargs) is different in g++ and VC when
you give arguments to run a program. This might be related to MPI as
well. The difference is one. Once the above changes are made, the
program is taking the correct argument.
However, it has been observed in the latest versions of sytem.cpp that
no modification needs be made to the file as distributed from the
LAMMPS website to work. The user however, instead of starting LAMMPS
by the command:
lammps in.file
as he would if he implemented the changes detailed here, would launch
in the Unix style:
lammps < in.file
7. The new version LAMMPS calls the error function:
double erfc(double)
This function is in the GNU C library. However, it's not found for
VC++.
Three options:
a. One can try to find erfc() from other libraries.
b. The erfc() is called for pair_modify table option. One can set
the table option to be 0 to avoid calling this function.
c. Write your own functions.
In this code, two files erfc.h, erfc.cpp are created and added to the project.
Files that call erfc() all add
#include "erfc.h" at the beginning.
Note: the functions are not fully tested, use with caution.
8. MSVC does not have a inttypes.h file. The simplest way
to deal with this problem is to download inttypes.h from the
following site:
http://www.koders.com/c/fidDE7D6EFFD475FAB1B7F6A2BBA791401CFA88FFA3.aspx
and add this file into the workspace.
9. MSVC does not have dirent.h. The problem is solved by downloading
a version of it for Windows from the following website:
http://www.softagalleria.net/dirent/index.en.html
10. Every time an error pops up for a line like this:
char *words[params_per_line];
replace to the following type:
char **words = new char*[params_per_line];
The dynamic memory allocation in MSVC requires a line like this.
11. Build the project. Specify appropriate input file to run the code.
The Windows result might be different from Unix results. Be Cautious.

1058
tools/ch2lmp/example/1ac7.pdb
View File

File diff suppressed because it is too large Load Diff

3744
tools/ch2lmp/example/1ac7.psf
View File

File diff suppressed because it is too large Load Diff

164
tools/lmp2cfg/README.txt
View File

@ -1,82 +1,82 @@
README.txt
lmp2cfg
Version 1.0
Ara Kooser
8/1/04
Contents
I. Before you start
II. Introduction
III. Compiling the program
IV. Running the program
I. Before you start
1) Read the READMEFIRST file
2) You will need either f77 or g77 compiler
3) A text editor (VI or emacs)
4) AtomEye( http://164.107.79.177/Archive/Graphics/A/) good as of 8/1/04
5) A sense of adventure and humor
II. Introduction
The lmp2cfg program converts a LAMMPS atom dump file into a series of .cfg
files.(000001.cfg, 000002.cfg,...). These files can be read into AtomEye
and made into a movie or viewed indivdually.
III. Compiling the program
The program should comply with either command f77 lmp2cfg.f or
g77 lmp2cfg.f
IV. Running the program
You need to create a user input file. In the examples folder you will find
an example input script, LAMMPS dump file, and the correct .cfg output.
The input script reads like this:
7 #number of atom types in your LAMMPS file
'dump.atom' #name of the LAMMPS dump file, you need the ' '
1 #first frame
10 #last frame
1 #first atom type
26.98154 #atomic weight
'ao' #atom name
2 #second atom type
15.9994 #atomic weight
'oh' #atom name
Make sure the input file and the atom dump file are in the same folder.
On the command line type
lmp2cfg < input_file.txt > out
You should get several .cfg files. For reading into AtomEye and making
movies see the AtomEye homepage.
If you get an error like:
open: 'new' file exists
apparent state: unit 21 named 00001.cfg
lately writing sequential formatted external IO
Abort
you need to first remove old copies of the cfg files.
If you get an error like:
open: No such file or directory
apparent state: unit 9 named dump.atom
last format: list io
lately reading direct formatted external IO
Abort
you need to check that the name of the dump file matches the name
in the input file, and that you enclosed it in single quotes.
README.txt
lmp2cfg
Version 1.0
Ara Kooser
8/1/04
Contents
I. Before you start
II. Introduction
III. Compiling the program
IV. Running the program
I. Before you start
1) Read the READMEFIRST file
2) You will need either f77 or g77 compiler
3) A text editor (VI or emacs)
4) AtomEye( http://164.107.79.177/Archive/Graphics/A/) good as of 8/1/04
5) A sense of adventure and humor
II. Introduction
The lmp2cfg program converts a LAMMPS atom dump file into a series of .cfg
files.(000001.cfg, 000002.cfg,...). These files can be read into AtomEye
and made into a movie or viewed indivdually.
III. Compiling the program
The program should comply with either command f77 lmp2cfg.f or
g77 lmp2cfg.f
IV. Running the program
You need to create a user input file. In the examples folder you will find
an example input script, LAMMPS dump file, and the correct .cfg output.
The input script reads like this:
7 #number of atom types in your LAMMPS file
'dump.atom' #name of the LAMMPS dump file, you need the ' '
1 #first frame
10 #last frame
1 #first atom type
26.98154 #atomic weight
'ao' #atom name
2 #second atom type
15.9994 #atomic weight
'oh' #atom name
Make sure the input file and the atom dump file are in the same folder.
On the command line type
lmp2cfg < input_file.txt > out
You should get several .cfg files. For reading into AtomEye and making
movies see the AtomEye homepage.
If you get an error like:
open: 'new' file exists
apparent state: unit 21 named 00001.cfg
lately writing sequential formatted external IO
Abort
you need to first remove old copies of the cfg files.
If you get an error like:
open: No such file or directory
apparent state: unit 9 named dump.atom
last format: list io
lately reading direct formatted external IO
Abort
you need to check that the name of the dump file matches the name
in the input file, and that you enclosed it in single quotes.

50
tools/lmp2cfg/examples/lmp2cfg_inp.txt
View File

@ -1,25 +1,25 @@
7
'dump.atom'
1
3
1
26.98154
'ao'
2
15.9994
'oh'
3
1.00794
'ho'
4
15.9994
'o*'
5
1.00794
'h*'
6
35.4527
'Cl'
7
207.2
'Pb'
7
'dump.atom'
1
3
1
26.98154
'ao'
2
15.9994
'oh'
3
1.00794
'ho'
4
15.9994
'o*'
5
1.00794
'h*'
6
35.4527
'Cl'
7
207.2
'Pb'

362
tools/lmp2cfg/lmp2cfg.f
View File

@ -1,181 +1,181 @@
*23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
c----------------------------------------------------------------------
program lammps to cfg
c----------------------------------------------------------------------
*Programming by Jeff Greathouse and Ara Kooser
*Version 1.0 9/1/04
*Sandia National Labs
*Converts LAMMPS atom dump to .cfg files for AtomEye
*This program is provided as is. Please see the README file.
c----------------------------------------------------------------------
implicit real*8 (a-h,o-z)
character*12 inhist,snapshot
character*4 fftype(100),name,ciframe
integer natom,iatom,itype(50000),q
integer atype(99),itycon,ntype,mass(99)
dimension x(50000),y(50000),z(50000)
dimension amass(99)
c-------Reads in the user input file-------------------------------
read(*,*) ntype
read(*,*) inhist
read(*,*) iframe1
read(*,*) iframe2
c write(*,*) ntype
c write(*,*) inhist
do 1, i=1, ntype
read(*,*) atype(i)
read(*,*) amass(i)
mass(i)=anint(amass(i))
read(*,*) fftype(i)
c write(*,*) atype(i)
c write(*,*) amass(i)
c write(*,*) fftype(i)
1 continue
c-------Lammps output file is 9, reads in lmps header--------------
name=inhist(1:4)
open(9,file=inhist,status='old',form='formatted')
c open(2,status='new',form='formatted')
iatom=0
iframe=0
jframe=0
c---------------------------------------------------------------------
c----------This begins the frame by frame reading section-------------
9999 continue
read(9,*,end=999)
read(9,*,end=999)
read(9,*,end=999)
read(9,*,end=999) natom
read(9,*,end=999)
read(9,50,end=999) xcell
read(9,50,end=999) ycell
read(9,50,end=999) zcell
read(9,*,end=999)
50 format(2x,f12.5)
1000 format(1x,i5,1x,i2,3f9.5)
do 440 j=1,natom
420 read(9,*,end=999)iatom,itype(iatom),x(iatom),y(iatom),z(iatom)
*23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
440 continue
jframe=jframe+1
if(jframe.lt.iframe1)goto 9999
iframe=iframe+1
if(iframe.ge.iframe2) goto 999
c--------------------------------------------------------------------
c-------This section writes each ts to a seperate .cfg file----------
c ciframe=char(iframe)
c snapshot(iframe)=name//ciframe//'.cfg'
c write(*,*)ciframe
c write(snapshot,'("Cfgs/",i7.7,".cfg")') iframe
ciframe='.cfg'
c write(*,*)ciframe
write(snapshot,'(i5.5,a4)')iframe,ciframe
c write(*,*)snapshot
open(unit=iframe+20,file=snapshot,status='new',
* form='formatted')
write((iframe+20),'(a22,i7)')'Number of particles = ',natom
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a16)')'A = 1.0 Angstrom'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),435)'H0(1,1) = ',xcell,' A'
write((iframe+20),'(a14)')'H0(1,2) = 0 A'
write((iframe+20),'(a14)')'H0(1,3) = 0 A'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a14)')'H0(2,1) = 0 A'
write((iframe+20),435)'H0(2,2) = ',ycell,' A'
write((iframe+20),'(a14)')'H0(2,3) = 0 A'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a14)')'H0(3,1) = 0 A'
write((iframe+20),'(a14)')'H0(3,2) = 0 A'
write((iframe+20),435)'H0(3,3) = ',zcell,' A'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a12)')'eta(1,1) = 0'
write((iframe+20),'(a12)')'eta(1,2) = 0'
write((iframe+20),'(a12)')'eta(1,3) = 0'
write((iframe+20),'(a12)')'eta(2,2) = 0'
write((iframe+20),'(a12)')'eta(2,3) = 0'
write((iframe+20),'(a12)')'eta(3,3) = 0'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
435 format(a11,f10.5,a2)
do 460, j=1,natom
do 450,i=1,ntype
c write(*,*)i,amass(i),fftype(i)
c write(*,*)ntype
if(itype(j).eq.atype(i))
* write((iframe+20),445)mass(i),fftype(i),x(j),
* y(j),z(j),' 0',' 0',' 0'
c---445 is the format for writing atom data to .cfg file------------
445 format(i3.3,1x,a2,1x,3f9.6,3a2)
450 continue
460 continue
go to 9999
999 continue
close(9)
end
*23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
c----------------------------------------------------------------------
program lammps to cfg
c----------------------------------------------------------------------
*Programming by Jeff Greathouse and Ara Kooser
*Version 1.0 9/1/04
*Sandia National Labs
*Converts LAMMPS atom dump to .cfg files for AtomEye
*This program is provided as is. Please see the README file.
c----------------------------------------------------------------------
implicit real*8 (a-h,o-z)
character*12 inhist,snapshot
character*4 fftype(100),name,ciframe
integer natom,iatom,itype(50000),q
integer atype(99),itycon,ntype,mass(99)
dimension x(50000),y(50000),z(50000)
dimension amass(99)
c-------Reads in the user input file-------------------------------
read(*,*) ntype
read(*,*) inhist
read(*,*) iframe1
read(*,*) iframe2
c write(*,*) ntype
c write(*,*) inhist
do 1, i=1, ntype
read(*,*) atype(i)
read(*,*) amass(i)
mass(i)=anint(amass(i))
read(*,*) fftype(i)
c write(*,*) atype(i)
c write(*,*) amass(i)
c write(*,*) fftype(i)
1 continue
c-------Lammps output file is 9, reads in lmps header--------------
name=inhist(1:4)
open(9,file=inhist,status='old',form='formatted')
c open(2,status='new',form='formatted')
iatom=0
iframe=0
jframe=0
c---------------------------------------------------------------------
c----------This begins the frame by frame reading section-------------
9999 continue
read(9,*,end=999)
read(9,*,end=999)
read(9,*,end=999)
read(9,*,end=999) natom
read(9,*,end=999)
read(9,50,end=999) xcell
read(9,50,end=999) ycell
read(9,50,end=999) zcell
read(9,*,end=999)
50 format(2x,f12.5)
1000 format(1x,i5,1x,i2,3f9.5)
do 440 j=1,natom
420 read(9,*,end=999)iatom,itype(iatom),x(iatom),y(iatom),z(iatom)
*23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
440 continue
jframe=jframe+1
if(jframe.lt.iframe1)goto 9999
iframe=iframe+1
if(iframe.ge.iframe2) goto 999
c--------------------------------------------------------------------
c-------This section writes each ts to a seperate .cfg file----------
c ciframe=char(iframe)
c snapshot(iframe)=name//ciframe//'.cfg'
c write(*,*)ciframe
c write(snapshot,'("Cfgs/",i7.7,".cfg")') iframe
ciframe='.cfg'
c write(*,*)ciframe
write(snapshot,'(i5.5,a4)')iframe,ciframe
c write(*,*)snapshot
open(unit=iframe+20,file=snapshot,status='new',
* form='formatted')
write((iframe+20),'(a22,i7)')'Number of particles = ',natom
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a16)')'A = 1.0 Angstrom'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),435)'H0(1,1) = ',xcell,' A'
write((iframe+20),'(a14)')'H0(1,2) = 0 A'
write((iframe+20),'(a14)')'H0(1,3) = 0 A'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a14)')'H0(2,1) = 0 A'
write((iframe+20),435)'H0(2,2) = ',ycell,' A'
write((iframe+20),'(a14)')'H0(2,3) = 0 A'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a14)')'H0(3,1) = 0 A'
write((iframe+20),'(a14)')'H0(3,2) = 0 A'
write((iframe+20),435)'H0(3,3) = ',zcell,' A'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a12)')'eta(1,1) = 0'
write((iframe+20),'(a12)')'eta(1,2) = 0'
write((iframe+20),'(a12)')'eta(1,3) = 0'
write((iframe+20),'(a12)')'eta(2,2) = 0'
write((iframe+20),'(a12)')'eta(2,3) = 0'
write((iframe+20),'(a12)')'eta(3,3) = 0'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),*)
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
write((iframe+20),'(a1)')'#'
435 format(a11,f10.5,a2)
do 460, j=1,natom
do 450,i=1,ntype
c write(*,*)i,amass(i),fftype(i)
c write(*,*)ntype
if(itype(j).eq.atype(i))
* write((iframe+20),445)mass(i),fftype(i),x(j),
* y(j),z(j),' 0',' 0',' 0'
c---445 is the format for writing atom data to .cfg file------------
445 format(i3.3,1x,a2,1x,3f9.6,3a2)
450 continue
460 continue
go to 9999
999 continue
close(9)
end

156
tools/lmp2traj/README.txt
View File

@ -1,78 +1,78 @@
README.txt
lmp2traj
Version 1.0
Ara Kooser
8/1/04
Contents
I. Before you start
II. Introduction
III. Compiling the program
IV. Running the program
I. Before you start
1) Read the READMEFIRST file
2) You will need either f77 or g77 compiler
3) A text editor (VI or emacs)
4) LAMMPS
5) A graphing program that can do contours
6) A sense of adventure and humor
II. Introduction
This program will take a LAMMPS atom dump file and provide the following three
files.
1) data for making contour maps
2) density profile
3) dipole information
III. Compiling the program
To compile the program run either ./f77 lmp2traj.f or
./g77 traj.f
IV. Running the program
First you need to create a user input file. There is an example of an input file
in the examples folder.
The input file reads like this:
'atom' # dump file name, needs the ' '
1 # first frame
60 # last frame
38.26119 # x dimension of the box
44.26119 # y dimension of the box
48.33150 # z dimension of the box
90. # angles of the box, always 90
90. # angles of the box, always 90
90. # angles of the box, always 90
82844.6 # volumne of the box in cubic Angstroms
5 # water oxygen atom type from LAMMPS (#)
6 # water hydrogen atom type from LAMMPS(#)
0. # leave at 0
5. # number of atom types
0. # z shift leave at 0
5. # number of density maps
'Surface (1) ho' # Enter name/description of atom
2 # atom type number from LAMMPS
'ho' # Column name for data
0 # Defines inner sphere, in A
48.3 # Defines out sphere, in A
Make sure you have the input file and the LAMMPS atom dump file in the same
directory.
To run the program type
lmp2traj < inp_file.txt > out
This should give you three files like in the examples/output folder.
README.txt
lmp2traj
Version 1.0
Ara Kooser
8/1/04
Contents
I. Before you start
II. Introduction
III. Compiling the program
IV. Running the program
I. Before you start
1) Read the READMEFIRST file
2) You will need either f77 or g77 compiler
3) A text editor (VI or emacs)
4) LAMMPS
5) A graphing program that can do contours
6) A sense of adventure and humor
II. Introduction
This program will take a LAMMPS atom dump file and provide the following three
files.
1) data for making contour maps
2) density profile
3) dipole information
III. Compiling the program
To compile the program run either ./f77 lmp2traj.f or
./g77 traj.f
IV. Running the program
First you need to create a user input file. There is an example of an input file
in the examples folder.
The input file reads like this:
'atom' # dump file name, needs the ' '
1 # first frame
60 # last frame
38.26119 # x dimension of the box
44.26119 # y dimension of the box
48.33150 # z dimension of the box
90. # angles of the box, always 90
90. # angles of the box, always 90
90. # angles of the box, always 90
82844.6 # volumne of the box in cubic Angstroms
5 # water oxygen atom type from LAMMPS (#)
6 # water hydrogen atom type from LAMMPS(#)
0. # leave at 0
5. # number of atom types
0. # z shift leave at 0
5. # number of density maps
'Surface (1) ho' # Enter name/description of atom
2 # atom type number from LAMMPS
'ho' # Column name for data
0 # Defines inner sphere, in A
48.3 # Defines out sphere, in A
Make sure you have the input file and the LAMMPS atom dump file in the same
directory.
To run the program type
lmp2traj < inp_file.txt > out
This should give you three files like in the examples/output folder.

52
tools/lmp2traj/examples/lmp2traj_inp.txt
View File

@ -1,26 +1,26 @@
'dump.atom'
1
3
29.00000
38.26120
44.79980
90.
90.
90.
49708.7291
1
2
0.
2.
0.
2
'Bulk (1) O-water'
1
'o*'
18.6
21.0
'Bulk (1) H-water'
2
'h*'
18.1
21.5
'dump.atom'
1
3
29.00000
38.26120
44.79980
90.
90.
90.
49708.7291
1
2
0.
2.
0.
2
'Bulk (1) O-water'
1
'o*'
18.6
21.0
'Bulk (1) H-water'
2
'h*'
18.1
21.5

504
tools/matlab/lmp2cfg.m
View File

@ -1,252 +1,252 @@
function lmp2cfg(varargin)
% Converts LAMMPS dump file to Extended CFG Format (No velocity) to be used
% with AtomEye (http://164.107.79.177/Archive/Graphics/A/)
% Input :
% Necessary (in order)
% timestep,Natoms,x_bound,y_bound,z_bound,H,atom_data,mass,
% cfg file name, dumpfile name
% Optional
% 'vel' --> 'yes' (if velocity data needs to be written)
% --> Default is 'no'
% 'atomtype' --> {Atom symbol}
% 'autonumber' --> 'yes' (default) numbers cfgfiles consecutively based on timestep
% 'autolog' --> 'on' (default) writes details of conversion to a log file
% 'rotation' --> To specify rotation (Transform in cfg file) say 'rotation' followed by
% number of rotations,the axes of rotation ('X' or 'Y' or 'Z')and the
% angle in degrees. CAUTION : Do remember that rotations are
% non-commutative and hence the order of rotations must be the same as
% intended
% 'aux' --> Auxiliary data must be a two column array with Quantity and Type in each
% column
%
%
% THE DATA MUST BE SCALED (0 to 1)IN ORDER FOR ATOMEYE TO READ PROPERLY
% Real coordinates x = s * H, x, s are 1x3 row vectors
%
% Example
% lmp2cfg_all(data.timestep,data.Natoms,data.xbound,data.ybound,...
% data.zbound,H,data.atom_data,mass,cfgfile,dumpfile,...
% 'vel','no','aux',{'sxx' 'Lammps O/P'
% 'syy' 'Lammps O/P'
% 'szz' 'Lammps O/P'
% 'sxy' 'Lammps O/P'
% 'sxz' 'Lammps O/P'
% 'syz' 'Lammps O/P'},...
% 'atomtype',{'Ni'
% 'Al'},...
% 'autonumber','on',...
% 'autolog','on',...
% 'rotation',2,'X',45,'Y',30);
%
% See also readdump_all, readdump_one, scandump
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
%------------ Defaults
vel = 'no';
auxstatus = 0;
atom_typestatus = 0;
rotstatus = 0;
autonumber_status = 1; % default is ON
autolog_status = 1;
%----------------------------------
%-----Required Input --------------
timestep = varargin{1}
Natoms = varargin{2};
x_bound = varargin{3};
y_bound = varargin{4};
z_bound = varargin{5};
H = varargin{6};
atom_data = varargin{7};
mass = varargin{8}; %arrange it in order of atom type , ie 1 - mass1 etc
filename = varargin{9}; %Just give one file name
if length(varargin) < 10
dumpfilename = ['dump.' filename];
else
dumpfilename = varargin{10};
end
if length(varargin) > 10
i=11;
while i<length(varargin)
id = varargin{i};
switch id
case 'vel'
vel = varargin{i+1};
i=i+2;
case 'aux'
auxiliary = varargin{i+1};
auxstatus = 1;
i=i+2;
case 'atomtype'
atom_type = varargin{i+1};
atom_typestatus = 1;
i=i+2;
case 'rotation'
nrot = varargin{i+1}; % number of rotations
i = i+2;
if nrot <=0
error('Number of rotations must be a positive value');
end
rotstatus = 1;
for j = 1 : 1 : nrot
rotaxes{j} = varargin{i};
i = i + 1;
rotangle(j) = varargin{i};
i = i + 1;
end
case 'autonumber'
numberstate = varargin{i+1};
i = i +2;
if strcmpi(numberstate,'ON')
autonumber_status = 1;
elseif strcmpi(numberstate,'OFF')
autonumber_status = 0;
end
case 'autolog'
logstate = varargin{i+1};
i = i + 2;
if strcmpi(logstate,'ON')
autolog_status = 1;
elseif strcmpi(logstate,'OFF')
autolog_status = 0;
end
end
end
end
% Calculating Transformation matrix [T]
if rotstatus == 0
T = eye(3); % Identity matrix
elseif rotstatus == 1
T = eye(3);
for j = 1 : 1 : nrot
B = beta(rotangle(j)*pi/180,rotaxes{j});
T = B*T;
end
end
T = T'; % because Transform = [beta]transpose
%----------------------------------
nfiles = length(timestep);
%----Default Atom type names------------
if atom_typestatus == 0
atom_type = {'Au'
'Ni'
'Zn'
'H'
'O'
'Cu'
'Al'
'Ag'
'C'
'Si'};
end
%--------Sorting atom types ---------
[s,id] = sort(atom_data(:,2,:));
%---------Writing CFG files---------------
for i = 1 : 1 : nfiles
if autonumber_status == 1
n = [filename '_' num2str(i) '.cfg'];
elseif autonumber_status == 0
n = [filename '.cfg'];
end
name{i}=n;
fid = fopen(name{i},'w+');
fprintf(fid,'Number of particles = %d\n',Natoms(i));
fprintf(fid,'A = 1.0000000000 Angstrom \n');
% Writing [H]
fprintf(fid,'H0(1,1) = %f A \n',H(1,1));
fprintf(fid,'H0(1,2) = %f A \n',H(1,2));
fprintf(fid,'H0(1,3) = %f A \n',H(1,3));
fprintf(fid,'H0(2,1) = %f A \n',H(2,1));
fprintf(fid,'H0(2,2) = %f A \n',H(2,2));
fprintf(fid,'H0(2,3) = %f A \n',H(2,3));
fprintf(fid,'H0(3,1) = %f A \n',H(3,1));
fprintf(fid,'H0(3,2) = %f A \n',H(3,2));
fprintf(fid,'H0(3,3) = %f A \n',H(3,3));
% Writing [T]
fprintf(fid,'Transform(1,1) = %f \n',T(1,1));
fprintf(fid,'Transform(1,2) = %f \n',T(1,2));
fprintf(fid,'Transform(1,3) = %f \n',T(1,3));
fprintf(fid,'Transform(2,1) = %f \n',T(2,1));
fprintf(fid,'Transform(2,2) = %f \n',T(2,2));
fprintf(fid,'Transform(2,3) = %f \n',T(2,3));
fprintf(fid,'Transform(3,1) = %f \n',T(3,1));
fprintf(fid,'Transform(3,2) = %f \n',T(3,2));
fprintf(fid,'Transform(3,3) = %f \n',T(3,3));
if strcmpi(vel,'no')
fprintf(fid,'.NO_VELOCITY. \n');
end
fprintf(fid,'entry_count = %d \n',length(atom_data(1,:,i))-2);
if auxstatus == 1
for k = 1 : 1 : length(auxiliary(:,1))
fprintf(fid,'auxiliary[%d] = %s [%s]\n',k-1,auxiliary{k,1},...
auxiliary{k,2});
end
end
aid = 1;
atom_change = 1;
for j = 1 : 1 : Natoms(i)
if atom_change == 1
fprintf(fid,'%f\n',mass(aid));
fprintf(fid,'%s\n',atom_type{aid});
end
atom_change = 0;
fprintf(fid,'%f\t',atom_data(id(j,1,i),3:length(atom_data(1,:,i)),i));
fprintf(fid,'\n');
if j ~= Natoms
if atom_data(id(j,1,i),2,i) ~= atom_data(id(j+1,1,i),2,i)
atom_change = 1;
aid = aid+1;
end
end
end
fclose(fid);
end
if autolog_status == 1
flog = fopen([filename '_lmp2cfg.log'],'w+');
fprintf(flog,'----------------------------------------------------------\n');
fprintf(flog,['LAMMPS DUMP to CFG file conversion :\t' datestr(now) '\n']);
fprintf(flog,'----------------------------------------------------------\n');
fprintf(flog,'LAMMPS Dump file : \t %s \n\n',dumpfilename);
for i = 1 : 1 : nfiles
fprintf(flog,'Timestep : %d --> \t\t %s \n',timestep(i),name{i});
end
fclose(flog);
end
%---------- Function to calculate beta
function b = beta(angle,axes)
switch axes
case 'X' % X axes
b = [1 0 0
0 cos(angle) sin(angle)
0 -sin(angle) cos(angle)];
case 'Y' % Y axes
b = [cos(angle) 0 sin(angle)
0 1 0
-sin(angle) 0 cos(angle)];
case 'Z' % Z axes
b = [cos(angle) sin(angle) 0
-sin(angle) cos(angle) 0
0 0 1];
end
end
% --------------------------------------------------
end % For main function
function lmp2cfg(varargin)
% Converts LAMMPS dump file to Extended CFG Format (No velocity) to be used
% with AtomEye (http://164.107.79.177/Archive/Graphics/A/)
% Input :
% Necessary (in order)
% timestep,Natoms,x_bound,y_bound,z_bound,H,atom_data,mass,
% cfg file name, dumpfile name
% Optional
% 'vel' --> 'yes' (if velocity data needs to be written)
% --> Default is 'no'
% 'atomtype' --> {Atom symbol}
% 'autonumber' --> 'yes' (default) numbers cfgfiles consecutively based on timestep
% 'autolog' --> 'on' (default) writes details of conversion to a log file
% 'rotation' --> To specify rotation (Transform in cfg file) say 'rotation' followed by
% number of rotations,the axes of rotation ('X' or 'Y' or 'Z')and the
% angle in degrees. CAUTION : Do remember that rotations are
% non-commutative and hence the order of rotations must be the same as
% intended
% 'aux' --> Auxiliary data must be a two column array with Quantity and Type in each
% column
%
%
% THE DATA MUST BE SCALED (0 to 1)IN ORDER FOR ATOMEYE TO READ PROPERLY
% Real coordinates x = s * H, x, s are 1x3 row vectors
%
% Example
% lmp2cfg_all(data.timestep,data.Natoms,data.xbound,data.ybound,...
% data.zbound,H,data.atom_data,mass,cfgfile,dumpfile,...
% 'vel','no','aux',{'sxx' 'Lammps O/P'
% 'syy' 'Lammps O/P'
% 'szz' 'Lammps O/P'
% 'sxy' 'Lammps O/P'
% 'sxz' 'Lammps O/P'
% 'syz' 'Lammps O/P'},...
% 'atomtype',{'Ni'
% 'Al'},...
% 'autonumber','on',...
% 'autolog','on',...
% 'rotation',2,'X',45,'Y',30);
%
% See also readdump_all, readdump_one, scandump
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
%------------ Defaults
vel = 'no';
auxstatus = 0;
atom_typestatus = 0;
rotstatus = 0;
autonumber_status = 1; % default is ON
autolog_status = 1;
%----------------------------------
%-----Required Input --------------
timestep = varargin{1}
Natoms = varargin{2};
x_bound = varargin{3};
y_bound = varargin{4};
z_bound = varargin{5};
H = varargin{6};
atom_data = varargin{7};
mass = varargin{8}; %arrange it in order of atom type , ie 1 - mass1 etc
filename = varargin{9}; %Just give one file name
if length(varargin) < 10
dumpfilename = ['dump.' filename];
else
dumpfilename = varargin{10};
end
if length(varargin) > 10
i=11;
while i<length(varargin)
id = varargin{i};
switch id
case 'vel'
vel = varargin{i+1};
i=i+2;
case 'aux'
auxiliary = varargin{i+1};
auxstatus = 1;
i=i+2;
case 'atomtype'
atom_type = varargin{i+1};
atom_typestatus = 1;
i=i+2;
case 'rotation'
nrot = varargin{i+1}; % number of rotations
i = i+2;
if nrot <=0
error('Number of rotations must be a positive value');
end
rotstatus = 1;
for j = 1 : 1 : nrot
rotaxes{j} = varargin{i};
i = i + 1;
rotangle(j) = varargin{i};
i = i + 1;
end
case 'autonumber'
numberstate = varargin{i+1};
i = i +2;
if strcmpi(numberstate,'ON')
autonumber_status = 1;
elseif strcmpi(numberstate,'OFF')
autonumber_status = 0;
end
case 'autolog'
logstate = varargin{i+1};
i = i + 2;
if strcmpi(logstate,'ON')
autolog_status = 1;
elseif strcmpi(logstate,'OFF')
autolog_status = 0;
end
end
end
end
% Calculating Transformation matrix [T]
if rotstatus == 0
T = eye(3); % Identity matrix
elseif rotstatus == 1
T = eye(3);
for j = 1 : 1 : nrot
B = beta(rotangle(j)*pi/180,rotaxes{j});
T = B*T;
end
end
T = T'; % because Transform = [beta]transpose
%----------------------------------
nfiles = length(timestep);
%----Default Atom type names------------
if atom_typestatus == 0
atom_type = {'Au'
'Ni'
'Zn'
'H'
'O'
'Cu'
'Al'
'Ag'
'C'
'Si'};
end
%--------Sorting atom types ---------
[s,id] = sort(atom_data(:,2,:));
%---------Writing CFG files---------------
for i = 1 : 1 : nfiles
if autonumber_status == 1
n = [filename '_' num2str(i) '.cfg'];
elseif autonumber_status == 0
n = [filename '.cfg'];
end
name{i}=n;
fid = fopen(name{i},'w+');
fprintf(fid,'Number of particles = %d\n',Natoms(i));
fprintf(fid,'A = 1.0000000000 Angstrom \n');
% Writing [H]
fprintf(fid,'H0(1,1) = %f A \n',H(1,1));
fprintf(fid,'H0(1,2) = %f A \n',H(1,2));
fprintf(fid,'H0(1,3) = %f A \n',H(1,3));
fprintf(fid,'H0(2,1) = %f A \n',H(2,1));
fprintf(fid,'H0(2,2) = %f A \n',H(2,2));
fprintf(fid,'H0(2,3) = %f A \n',H(2,3));
fprintf(fid,'H0(3,1) = %f A \n',H(3,1));
fprintf(fid,'H0(3,2) = %f A \n',H(3,2));
fprintf(fid,'H0(3,3) = %f A \n',H(3,3));
% Writing [T]
fprintf(fid,'Transform(1,1) = %f \n',T(1,1));
fprintf(fid,'Transform(1,2) = %f \n',T(1,2));
fprintf(fid,'Transform(1,3) = %f \n',T(1,3));
fprintf(fid,'Transform(2,1) = %f \n',T(2,1));
fprintf(fid,'Transform(2,2) = %f \n',T(2,2));
fprintf(fid,'Transform(2,3) = %f \n',T(2,3));
fprintf(fid,'Transform(3,1) = %f \n',T(3,1));
fprintf(fid,'Transform(3,2) = %f \n',T(3,2));
fprintf(fid,'Transform(3,3) = %f \n',T(3,3));
if strcmpi(vel,'no')
fprintf(fid,'.NO_VELOCITY. \n');
end
fprintf(fid,'entry_count = %d \n',length(atom_data(1,:,i))-2);
if auxstatus == 1
for k = 1 : 1 : length(auxiliary(:,1))
fprintf(fid,'auxiliary[%d] = %s [%s]\n',k-1,auxiliary{k,1},...
auxiliary{k,2});
end
end
aid = 1;
atom_change = 1;
for j = 1 : 1 : Natoms(i)
if atom_change == 1
fprintf(fid,'%f\n',mass(aid));
fprintf(fid,'%s\n',atom_type{aid});
end
atom_change = 0;
fprintf(fid,'%f\t',atom_data(id(j,1,i),3:length(atom_data(1,:,i)),i));
fprintf(fid,'\n');
if j ~= Natoms
if atom_data(id(j,1,i),2,i) ~= atom_data(id(j+1,1,i),2,i)
atom_change = 1;
aid = aid+1;
end
end
end
fclose(fid);
end
if autolog_status == 1
flog = fopen([filename '_lmp2cfg.log'],'w+');
fprintf(flog,'----------------------------------------------------------\n');
fprintf(flog,['LAMMPS DUMP to CFG file conversion :\t' datestr(now) '\n']);
fprintf(flog,'----------------------------------------------------------\n');
fprintf(flog,'LAMMPS Dump file : \t %s \n\n',dumpfilename);
for i = 1 : 1 : nfiles
fprintf(flog,'Timestep : %d --> \t\t %s \n',timestep(i),name{i});
end
fclose(flog);
end
%---------- Function to calculate beta
function b = beta(angle,axes)
switch axes
case 'X' % X axes
b = [1 0 0
0 cos(angle) sin(angle)
0 -sin(angle) cos(angle)];
case 'Y' % Y axes
b = [cos(angle) 0 sin(angle)
0 1 0
-sin(angle) 0 cos(angle)];
case 'Z' % Z axes
b = [cos(angle) sin(angle) 0
-sin(angle) cos(angle) 0
0 0 1];
end
end
% --------------------------------------------------
end % For main function

332
tools/matlab/readEAM.m
View File

@ -1,166 +1,166 @@
function varargout = readEAM(varargin)
% Function to read EAM potential files
% Input
% 1: EAM potential file name
% 2: file type --> 'FUNCFL' for single element file
% --> 'SETFL' for multiple element file
% Output : Structure with members
% .embed : Embedding function
% .pair : Pair Potential
% .elecden : Electron Density
% .nrho : Number of points for electron density
% .drho : increment of r for electron density
% .nr : Number of points for pair potential
% .dr : increment of r for pair potential
% .rcut : cut-off distance
% All output is in exactly the same units as in the EAM file. For
% multielement SETFL files the members are multidimensional arrays with
% each element data stored in (:,:,ELEM)
%
% Example
% eam = readEAM('cuu3.eam','funcfl');
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
if length(varargin) < 2
error('Too few input parameters : Need Filename and filetype');
elseif length(varargin) > 2
error('Too many input parameters');
end
filename = varargin{1};
type = varargin{2};
try
fid = fopen(filename,'r');
catch
error('EAM file not found!');
end
if strcmpi(type,'FUNCFL')
t = 1;
elseif strcmpi(type,'SETFL')
t = 2;
else
error(['Unknown file type : "' type '"']);
end
switch t
case 1 %FUNCFL
% Read line 1 : comment line
a = fgetl(fid);
% read Line 2
rem = fgetl(fid); % ielem amass blat lat
[ielem,rem] = strtok(rem);
[amass,rem] = strtok(rem);
[blat,rem] = strtok(rem);
[lat,rem] = strtok(rem);
ielem = str2num(ielem); % Atomic Number
amass = str2num(amass); % Atomic Mass
blat = str2num(blat); % Lattice Constant
% Read Line 3
a = str2num(fgetl(fid));
nrho = a(1);
drho = a(2);
nr = a(3);
dr = a(4);
rcut = a(5);
% Reading embedding function
for i = 1 : 1 : nrho/5
embed(i,:) = str2num(fgetl(fid));
end
% Reading pair potential
for i = 1 : 1 : nr/5
pair(i,:) = str2num(fgetl(fid));
end
% Reading electron density function
for i = 1 : 1 : nr/5
elecden(i,:) = str2num(fgetl(fid));
end
% Output
out.embed = embed;
out.pair = pair;
out.elecden = elecden;
out.nrho = nrho;
out.drho = drho;
out.nr = nr;
out.dr = dr;
out.rcut = rcut;
varargout{1} = out;
% --------------------------------------------------------------------
case 2 % SETFL
% Read lines 1 - 3 : comment lines
a = fgetl(fid);
a = fgetl(fid);
a = fgetl(fid);
% Atom types
ntypes = str2num(fgetl(fid));
% Read Global information
a = str2num(fgetl(fid));
nrho = a(1);
drho = a(2);
nr = a(3);
dr = a(4);
rcut = a(5);
% Read element specific Data
% Embedding function and Electron Density
for elem = 1 : 1 : ntypes
rem = fgetl(fid); % ielem amass blat lat
[ielem1,rem] = strtok(rem);
[amass1,rem] = strtok(rem);
[blat1,rem] = strtok(rem);
[lat1,rem] = strtok(rem);
ielem(elem) = str2num(ielem1); % Atomic Number
amass(elem) = str2num(amass1); % Atomic Mass
blat(elem) = str2num(blat1); % Lattice Constant
lat(elem,:) = lat1; % Lattice type
% Reading embedding function
for i = 1 : 1 : nrho/5
embed(i,:,elem) = str2num(fgetl(fid));
end
% Reading electron density function
for i = 1 : 1 : nr/5
elecden(i,:,elem) = str2num(fgetl(fid));
end
end
% Pair Potentials
n_pair = ntypes + (factorial(ntypes)/2);
for np = 1 : 1 : n_pair
for i = 1 : 1 : nr/5
pair(i,:,np) = str2num(fgetl(fid));
end
end
% Output
out.embed = embed;
out.elecden = elecden;
out.pair = pair;
out.nrho = nrho;
out.drho = drho;
out.nr = nr;
out.dr = dr;
out.rcut = rcut;
varargout{1} = out;
end
function varargout = readEAM(varargin)
% Function to read EAM potential files
% Input
% 1: EAM potential file name
% 2: file type --> 'FUNCFL' for single element file
% --> 'SETFL' for multiple element file
% Output : Structure with members
% .embed : Embedding function
% .pair : Pair Potential
% .elecden : Electron Density
% .nrho : Number of points for electron density
% .drho : increment of r for electron density
% .nr : Number of points for pair potential
% .dr : increment of r for pair potential
% .rcut : cut-off distance
% All output is in exactly the same units as in the EAM file. For
% multielement SETFL files the members are multidimensional arrays with
% each element data stored in (:,:,ELEM)
%
% Example
% eam = readEAM('cuu3.eam','funcfl');
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
if length(varargin) < 2
error('Too few input parameters : Need Filename and filetype');
elseif length(varargin) > 2
error('Too many input parameters');
end
filename = varargin{1};
type = varargin{2};
try
fid = fopen(filename,'r');
catch
error('EAM file not found!');
end
if strcmpi(type,'FUNCFL')
t = 1;
elseif strcmpi(type,'SETFL')
t = 2;
else
error(['Unknown file type : "' type '"']);
end
switch t
case 1 %FUNCFL
% Read line 1 : comment line
a = fgetl(fid);
% read Line 2
rem = fgetl(fid); % ielem amass blat lat
[ielem,rem] = strtok(rem);
[amass,rem] = strtok(rem);
[blat,rem] = strtok(rem);
[lat,rem] = strtok(rem);
ielem = str2num(ielem); % Atomic Number
amass = str2num(amass); % Atomic Mass
blat = str2num(blat); % Lattice Constant
% Read Line 3
a = str2num(fgetl(fid));
nrho = a(1);
drho = a(2);
nr = a(3);
dr = a(4);
rcut = a(5);
% Reading embedding function
for i = 1 : 1 : nrho/5
embed(i,:) = str2num(fgetl(fid));
end
% Reading pair potential
for i = 1 : 1 : nr/5
pair(i,:) = str2num(fgetl(fid));
end
% Reading electron density function
for i = 1 : 1 : nr/5
elecden(i,:) = str2num(fgetl(fid));
end
% Output
out.embed = embed;
out.pair = pair;
out.elecden = elecden;
out.nrho = nrho;
out.drho = drho;
out.nr = nr;
out.dr = dr;
out.rcut = rcut;
varargout{1} = out;
% --------------------------------------------------------------------
case 2 % SETFL
% Read lines 1 - 3 : comment lines
a = fgetl(fid);
a = fgetl(fid);
a = fgetl(fid);
% Atom types
ntypes = str2num(fgetl(fid));
% Read Global information
a = str2num(fgetl(fid));
nrho = a(1);
drho = a(2);
nr = a(3);
dr = a(4);
rcut = a(5);
% Read element specific Data
% Embedding function and Electron Density
for elem = 1 : 1 : ntypes
rem = fgetl(fid); % ielem amass blat lat
[ielem1,rem] = strtok(rem);
[amass1,rem] = strtok(rem);
[blat1,rem] = strtok(rem);
[lat1,rem] = strtok(rem);
ielem(elem) = str2num(ielem1); % Atomic Number
amass(elem) = str2num(amass1); % Atomic Mass
blat(elem) = str2num(blat1); % Lattice Constant
lat(elem,:) = lat1; % Lattice type
% Reading embedding function
for i = 1 : 1 : nrho/5
embed(i,:,elem) = str2num(fgetl(fid));
end
% Reading electron density function
for i = 1 : 1 : nr/5
elecden(i,:,elem) = str2num(fgetl(fid));
end
end
% Pair Potentials
n_pair = ntypes + (factorial(ntypes)/2);
for np = 1 : 1 : n_pair
for i = 1 : 1 : nr/5
pair(i,:,np) = str2num(fgetl(fid));
end
end
% Output
out.embed = embed;
out.elecden = elecden;
out.pair = pair;
out.nrho = nrho;
out.drho = drho;
out.nr = nr;
out.dr = dr;
out.rcut = rcut;
varargout{1} = out;
end

240
tools/matlab/readdump_one.m
View File

@ -1,120 +1,120 @@
function [varargout] = readdump_one(varargin)
% Read LAMMPS dump file one timestep at a time
% Input
% Dump file name with path
% Starting file pointer position in dump file
% Number of columns in the dump file
% Output is in the form of a structure with following variables
% .timestep --> Vector containing all time steps
% .Natoms --> Vector containing number of atoms at each time step
% .x_bound --> [t,2] array with xlo,xhi at each time step
% .y_bound --> [t,2] array with ylo,yhi at each time step
% .z_bound --> [t,2] array with zlo,zhi at each time step
% .atom_data --> 2 dimensional array with data
% .position --> file pointer for reading next time
%
% Example
% data = readdump_one('dump.LAMMPS',0,5);
% Reads the first timestep in the file dump.LAMMPS
% Specify position = -1 if only the last dump step in the
% file is needed
%
% See also readdump, scandump
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
try
dump = fopen(varargin{1},'r');
catch
error('Dumpfile not found!');
end
position = varargin{2}; % from beg of file
ncol = varargin{3}; %number of columns
i=1;
t=0;
done = 0;
last_status = 0;
if position ~= -1
fseek(dump,position,'bof');
else
last_status = 1;
end
while done == 0 & last_status == 0
id = fgetl(dump);
switch id
case 'ITEM: TIMESTEP'
if t == 0
timestep(i) = str2num(fgetl(dump));
t=1;
end
case 'ITEM: NUMBER OF ATOMS'
Natoms = str2num(fgetl(dump));
case 'ITEM: BOX BOUNDS'
x_bound(1,:) = str2num(fgetl(dump));
y_bound(1,:) = str2num(fgetl(dump));
z_bound(1,:) = str2num(fgetl(dump));
case 'ITEM: ATOMS'
atom_data = zeros(Natoms,ncol);%Allocate memory for atom data
for j = 1 : 1: Natoms
atom_data(j,:) = str2num(fgetl(dump));
end
done = 1;
p = ftell(dump);
end
end
% Getting only the last step
if last_status == 1
% First get the position of the beginning of the last step in the
% dumpfile
while ~feof(dump)
temp = fgetl(dump);
if length(temp) == 14
if strcmpi(temp,'ITEM: TIMESTEP')
p = ftell(dump); % starting position of line next to the header
end
end
end
fclose(dump);
dump = fopen(varargin{1},'r');
fseek(dump,p,'bof');
% getting Timestep
timestep = str2num(fgetl(dump));
while ~feof(dump)
id = fgetl(dump);
switch id
case 'ITEM: NUMBER OF ATOMS'
Natoms = str2num(fgetl(dump));
case 'ITEM: BOX BOUNDS'
x_bound(1,:) = str2num(fgetl(dump));
y_bound(1,:) = str2num(fgetl(dump));
z_bound(1,:) = str2num(fgetl(dump));
case 'ITEM: ATOMS'
atom_data = zeros(Natoms,ncol);%Allocate memory for atom data
for j = 1 : 1: Natoms
atom_data(j,:) = str2num(fgetl(dump));
end
end
end
end
%----------Outputs-------------
%OUTPUTS IN SAME VARIABLE STRUCTURE
varargout{1}.timestep = timestep;
varargout{1}.Natoms = Natoms;
varargout{1}.x_bound = x_bound;
varargout{1}.y_bound = y_bound;
varargout{1}.z_bound = z_bound;
varargout{1}.atom_data = atom_data;
varargout{1}.position = p; %gives postion of ITEM: TIMESTEP line
%------------------------------
fclose(dump);
function [varargout] = readdump_one(varargin)
% Read LAMMPS dump file one timestep at a time
% Input
% Dump file name with path
% Starting file pointer position in dump file
% Number of columns in the dump file
% Output is in the form of a structure with following variables
% .timestep --> Vector containing all time steps
% .Natoms --> Vector containing number of atoms at each time step
% .x_bound --> [t,2] array with xlo,xhi at each time step
% .y_bound --> [t,2] array with ylo,yhi at each time step
% .z_bound --> [t,2] array with zlo,zhi at each time step
% .atom_data --> 2 dimensional array with data
% .position --> file pointer for reading next time
%
% Example
% data = readdump_one('dump.LAMMPS',0,5);
% Reads the first timestep in the file dump.LAMMPS
% Specify position = -1 if only the last dump step in the
% file is needed
%
% See also readdump, scandump
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
try
dump = fopen(varargin{1},'r');
catch
error('Dumpfile not found!');
end
position = varargin{2}; % from beg of file
ncol = varargin{3}; %number of columns
i=1;
t=0;
done = 0;
last_status = 0;
if position ~= -1
fseek(dump,position,'bof');
else
last_status = 1;
end
while done == 0 & last_status == 0
id = fgetl(dump);
switch id
case 'ITEM: TIMESTEP'
if t == 0
timestep(i) = str2num(fgetl(dump));
t=1;
end
case 'ITEM: NUMBER OF ATOMS'
Natoms = str2num(fgetl(dump));
case 'ITEM: BOX BOUNDS'
x_bound(1,:) = str2num(fgetl(dump));
y_bound(1,:) = str2num(fgetl(dump));
z_bound(1,:) = str2num(fgetl(dump));
case 'ITEM: ATOMS'
atom_data = zeros(Natoms,ncol);%Allocate memory for atom data
for j = 1 : 1: Natoms
atom_data(j,:) = str2num(fgetl(dump));
end
done = 1;
p = ftell(dump);
end
end
% Getting only the last step
if last_status == 1
% First get the position of the beginning of the last step in the
% dumpfile
while ~feof(dump)
temp = fgetl(dump);
if length(temp) == 14
if strcmpi(temp,'ITEM: TIMESTEP')
p = ftell(dump); % starting position of line next to the header
end
end
end
fclose(dump);
dump = fopen(varargin{1},'r');
fseek(dump,p,'bof');
% getting Timestep
timestep = str2num(fgetl(dump));
while ~feof(dump)
id = fgetl(dump);
switch id
case 'ITEM: NUMBER OF ATOMS'
Natoms = str2num(fgetl(dump));
case 'ITEM: BOX BOUNDS'
x_bound(1,:) = str2num(fgetl(dump));
y_bound(1,:) = str2num(fgetl(dump));
z_bound(1,:) = str2num(fgetl(dump));
case 'ITEM: ATOMS'
atom_data = zeros(Natoms,ncol);%Allocate memory for atom data
for j = 1 : 1: Natoms
atom_data(j,:) = str2num(fgetl(dump));
end
end
end
end
%----------Outputs-------------
%OUTPUTS IN SAME VARIABLE STRUCTURE
varargout{1}.timestep = timestep;
varargout{1}.Natoms = Natoms;
varargout{1}.x_bound = x_bound;
varargout{1}.y_bound = y_bound;
varargout{1}.z_bound = z_bound;
varargout{1}.atom_data = atom_data;
varargout{1}.position = p; %gives postion of ITEM: TIMESTEP line
%------------------------------
fclose(dump);

182
tools/matlab/readlog.m
View File

@ -1,92 +1,92 @@
function [varargout] = readlog(varargin)
% Read LAMMPS log files
% input is log file name with path
% output is a structure --> out
% out.Chead --> has heading of columns in each run
% out.data --> has the data during each run in the form of string
% Type str2num(out.data{i}) to get the numeric array
%
% Example
% logdata = readlog('log.LAMMPS');
%
% Author : sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% Arun K. Subramaniyan
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
logfile = varargin{1};
try
fid = fopen(logfile,'r');
catch
error('Log file not found!');
end
loop = 1;
while feof(fid) == 0
%----------- To get first line of thermo output --------
while feof(fid) == 0
a = fgetl(fid);
if length(a) > 4 && strcmp(a(1:4),'Step')
findstep = 1;
break;
end
end
%-------------------------------------------------------
%---------Seperate column headings----------------------
j=1;
k=1;
Ch='';
i=1;
while i < length(a) && findstep == 1
for i = k : 1 : length(a)
if strcmp(a(i),' ')
Chead{loop,j} = Ch;
clear Ch;
Ch = '';
k=i+1;
j=j+1;
break;
else
Ch = [Ch a(i)];
end
end
end
if findstep == 1
Chead{loop,j} = Ch; % to get the last column heading
end
%-------------------------------------------------------
%----------------------Get Data-------------------------
id = 1; %row id...
while feof(fid) == 0
a = fgetl(fid);
if strcmp(a(1:4),'Loop')
loop = loop + 1;
break;
else
logdata(id,:) = str2num(a);
id = id+1;
end
end
%--------------------------------------------------------
if feof(fid) == 0
b = num2str(logdata);
data{loop-1} = b;
clear b;
clear logdata;
findstep = 0;
end
end
fclose(fid);
%--------OUTPUT-------------------------------------------
out.Chead = Chead;
out.data = data;
varargout{1} = out;
function [varargout] = readlog(varargin)
% Read LAMMPS log files
% input is log file name with path
% output is a structure --> out
% out.Chead --> has heading of columns in each run
% out.data --> has the data during each run in the form of string
% Type str2num(out.data{i}) to get the numeric array
%
% Example
% logdata = readlog('log.LAMMPS');
%
% Author : sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% Arun K. Subramaniyan
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
logfile = varargin{1};
try
fid = fopen(logfile,'r');
catch
error('Log file not found!');
end
loop = 1;
while feof(fid) == 0
%----------- To get first line of thermo output --------
while feof(fid) == 0
a = fgetl(fid);
if length(a) > 4 && strcmp(a(1:4),'Step')
findstep = 1;
break;
end
end
%-------------------------------------------------------
%---------Seperate column headings----------------------
j=1;
k=1;
Ch='';
i=1;
while i < length(a) && findstep == 1
for i = k : 1 : length(a)
if strcmp(a(i),' ')
Chead{loop,j} = Ch;
clear Ch;
Ch = '';
k=i+1;
j=j+1;
break;
else
Ch = [Ch a(i)];
end
end
end
if findstep == 1
Chead{loop,j} = Ch; % to get the last column heading
end
%-------------------------------------------------------
%----------------------Get Data-------------------------
id = 1; %row id...
while feof(fid) == 0
a = fgetl(fid);
if strcmp(a(1:4),'Loop')
loop = loop + 1;
break;
else
logdata(id,:) = str2num(a);
id = id+1;
end
end
%--------------------------------------------------------
if feof(fid) == 0
b = num2str(logdata);
data{loop-1} = b;
clear b;
clear logdata;
findstep = 0;
end
end
fclose(fid);
%--------OUTPUT-------------------------------------------
out.Chead = Chead;
out.data = data;
varargout{1} = out;

284
tools/matlab/readrdf.m
View File

@ -1,142 +1,142 @@
function varargout = readrdf(varargin)
% Function to read Radial Distribution Funtion output from LAMMPS
% Input
% 'bin' --> number of bins in rdf histogram
% 'runtime' --> Run length of each of the run commands
% 'step' --> rdf Dump step for each of the run commands
% 'ncol' --> number of columns in the file
% 'final' --> 'YES' indicates that only the average values will be extracted
% If only the averages are needed, you don't need to specify 'runtime',
% 'step' and 'ncol'
%
% Output is in the form of a structure with following variables
%
% if 'final' was given as 'NO'
% .rdf_step_data --> Matrix with RDF for each of the dumped steps
% .rdf_ave_data --> Matrix with average RDF for each RUN
% .rdf_ave_data --> Matrix with average RDF for each RUN
%
% if 'final' was given as 'YES'
% .rdf_ave_data --> Matrix with average RDF for each RUN
%
% Example
% rdf = readrdf('inputfile','bin',100,'runtime',[30000;20000],...
% 'step',[100;100],'ncol',3,'final','no')
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
rdf_name = varargin{1}; % RDF File name
% Setting all status checks to zero
bin_status = 0;
runtime_status = 0;
step_status = 0;
ncol_status = 0;
final_status = 'no';
for id = 2 : 1 : length(varargin)
if strcmpi(varargin{id},'bin')
rdf_bin = varargin{id+1}; % RDF bins
bin_status = 1;
elseif strcmpi(varargin{id},'runtime')
rdf_runtime = varargin{id+1}; % Runtimes
runtime_status = 1;
elseif strcmpi(varargin{id},'step')
rdf_step = varargin{id+1}; % Runtimes
step_status = 1;
elseif strcmpi(varargin{id},'ncol')
ncol = varargin{id+1}; % Runtimes
ncol_status = 1;
elseif strcmpi(varargin{id},'final')
final_status = varargin{id+1};
end
end
if ncol_status == 0
ncol = 3;
end
% Check for errors in input arguments
if bin_status == 0
error('No bin specified');
elseif step_status == 1 && runtime_status == 0
error('Step size specified without Runtime');
elseif step_status == 0 && runtime_status == 1
error('Runtime specified without Step size');
end
if step_status == 1 && runtime_status == 1
if length(rdf_runtime) ~= length(rdf_step)
error('Runtime and Step size do not match');
end
end
% Preallocating memory if runtime and step size are provided
if step_status == 1 && runtime_status == 1 && strcmpi(final_status,'no')
total_steps = round(sum(rdf_runtime./rdf_step));
rdf_step_data = zeros(rdf_bin,ncol,total_steps);
rdf_ave_data = zeros(rdf_bin,ncol,length(rdf_runtime));
elseif strcmpi(final_status,'yes') && step_status == 1 && runtime_status == 1
rdf_ave_data = zeros(rdf_bin,ncol,length(rdf_runtime));
end
try
rdf_file = fopen(rdf_name,'r');
catch
error('RDF file not found!');
end
if strcmpi(final_status,'yes')
run_id = 1; % Run id..
while feof(rdf_file) ~= 1
rdf_data = fgetl(rdf_file);
if strcmpi(rdf_data(1:11),'RUN AVERAGE')
fgetl(rdf_file); % to skip the title
for id = 1 : 1 : rdf_bin
rdf_ave_data(id,:,run_id) = str2num(fgetl(rdf_file));
end
run_id = run_id + 1;
end
end
else
run_id = 1;
id = 1;
while feof(rdf_file) ~= 1
rdf_data = fgetl(rdf_file);
if strcmpi(rdf_data(1:8),'TIMESTEP')
timestep(id,1) = str2num(rdf_data(10:length(rdf_data)));
fgetl(rdf_file); % to skip the title
for j = 1 : 1 : rdf_bin
rdf_step_data(j,:,id) = str2num(fgetl(rdf_file));
end
id = id+1;
elseif strcmpi(rdf_data(1:11),'RUN AVERAGE')
fgetl(rdf_file); % to skip the title
for j = 1 : 1 : rdf_bin
rdf_ave_data(j,:,run_id) = str2num(fgetl(rdf_file));
end
run_id = run_id + 1;
end
end
end
fclose(rdf_file);
if strcmpi(final_status,'no')
out_data.timestep = timestep;
out_data.rdf_step_data = rdf_step_data;
out_data.rdf_ave_data = rdf_ave_data;
else
out_data.rdf_ave_data = rdf_ave_data;
end
varargout{1} = out_data;
function varargout = readrdf(varargin)
% Function to read Radial Distribution Funtion output from LAMMPS
% Input
% 'bin' --> number of bins in rdf histogram
% 'runtime' --> Run length of each of the run commands
% 'step' --> rdf Dump step for each of the run commands
% 'ncol' --> number of columns in the file
% 'final' --> 'YES' indicates that only the average values will be extracted
% If only the averages are needed, you don't need to specify 'runtime',
% 'step' and 'ncol'
%
% Output is in the form of a structure with following variables
%
% if 'final' was given as 'NO'
% .rdf_step_data --> Matrix with RDF for each of the dumped steps
% .rdf_ave_data --> Matrix with average RDF for each RUN
% .rdf_ave_data --> Matrix with average RDF for each RUN
%
% if 'final' was given as 'YES'
% .rdf_ave_data --> Matrix with average RDF for each RUN
%
% Example
% rdf = readrdf('inputfile','bin',100,'runtime',[30000;20000],...
% 'step',[100;100],'ncol',3,'final','no')
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
rdf_name = varargin{1}; % RDF File name
% Setting all status checks to zero
bin_status = 0;
runtime_status = 0;
step_status = 0;
ncol_status = 0;
final_status = 'no';
for id = 2 : 1 : length(varargin)
if strcmpi(varargin{id},'bin')
rdf_bin = varargin{id+1}; % RDF bins
bin_status = 1;
elseif strcmpi(varargin{id},'runtime')
rdf_runtime = varargin{id+1}; % Runtimes
runtime_status = 1;
elseif strcmpi(varargin{id},'step')
rdf_step = varargin{id+1}; % Runtimes
step_status = 1;
elseif strcmpi(varargin{id},'ncol')
ncol = varargin{id+1}; % Runtimes
ncol_status = 1;
elseif strcmpi(varargin{id},'final')
final_status = varargin{id+1};
end
end
if ncol_status == 0
ncol = 3;
end
% Check for errors in input arguments
if bin_status == 0
error('No bin specified');
elseif step_status == 1 && runtime_status == 0
error('Step size specified without Runtime');
elseif step_status == 0 && runtime_status == 1
error('Runtime specified without Step size');
end
if step_status == 1 && runtime_status == 1
if length(rdf_runtime) ~= length(rdf_step)
error('Runtime and Step size do not match');
end
end
% Preallocating memory if runtime and step size are provided
if step_status == 1 && runtime_status == 1 && strcmpi(final_status,'no')
total_steps = round(sum(rdf_runtime./rdf_step));
rdf_step_data = zeros(rdf_bin,ncol,total_steps);
rdf_ave_data = zeros(rdf_bin,ncol,length(rdf_runtime));
elseif strcmpi(final_status,'yes') && step_status == 1 && runtime_status == 1
rdf_ave_data = zeros(rdf_bin,ncol,length(rdf_runtime));
end
try
rdf_file = fopen(rdf_name,'r');
catch
error('RDF file not found!');
end
if strcmpi(final_status,'yes')
run_id = 1; % Run id..
while feof(rdf_file) ~= 1
rdf_data = fgetl(rdf_file);
if strcmpi(rdf_data(1:11),'RUN AVERAGE')
fgetl(rdf_file); % to skip the title
for id = 1 : 1 : rdf_bin
rdf_ave_data(id,:,run_id) = str2num(fgetl(rdf_file));
end
run_id = run_id + 1;
end
end
else
run_id = 1;
id = 1;
while feof(rdf_file) ~= 1
rdf_data = fgetl(rdf_file);
if strcmpi(rdf_data(1:8),'TIMESTEP')
timestep(id,1) = str2num(rdf_data(10:length(rdf_data)));
fgetl(rdf_file); % to skip the title
for j = 1 : 1 : rdf_bin
rdf_step_data(j,:,id) = str2num(fgetl(rdf_file));
end
id = id+1;
elseif strcmpi(rdf_data(1:11),'RUN AVERAGE')
fgetl(rdf_file); % to skip the title
for j = 1 : 1 : rdf_bin
rdf_ave_data(j,:,run_id) = str2num(fgetl(rdf_file));
end
run_id = run_id + 1;
end
end
end
fclose(rdf_file);
if strcmpi(final_status,'no')
out_data.timestep = timestep;
out_data.rdf_step_data = rdf_step_data;
out_data.rdf_ave_data = rdf_ave_data;
else
out_data.rdf_ave_data = rdf_ave_data;
end
varargout{1} = out_data;

6
tools/msi2lmp/src/GetParameters.c
View File

@ -261,7 +261,7 @@ void GetParameters(int Forcefield)
for (j=0; j < 3; j++)
fprintf(stderr," %-3s", atomtypes[angletypes[i].types[j]].potential);
for (j=0; j < 4; j++) fprintf(stderr," %8.4f",angletypes[i].params[j]);
fprintf(stderr,"\n");
fprintf(stderr,"\n");
}
if (forcefield > 1) {
@ -448,7 +448,7 @@ void GetParameters(int Forcefield)
}
else {
dihedraltypes[i].midbonddihedral_cross_term[0] =
ff_midbontor.data[k].ff_param[0];
ff_midbontor.data[k].ff_param[0];
dihedraltypes[i].midbonddihedral_cross_term[1] =
ff_midbontor.data[k].ff_param[1];
dihedraltypes[i].midbonddihedral_cross_term[2] =
@ -482,7 +482,7 @@ void GetParameters(int Forcefield)
}
else {
dihedraltypes[i].angledihedral_cross_term[0] =
ff_angtor.data[k].ff_param[0];
ff_angtor.data[k].ff_param[0];
dihedraltypes[i].angledihedral_cross_term[1] =
ff_angtor.data[k].ff_param[1];
dihedraltypes[i].angledihedral_cross_term[2] =

2
tools/msi2lmp/src/InitializeItems.c
View File

@ -33,7 +33,7 @@ void InitializeItems(void)
}
else {
strcpy(ff_bond.keyword, "#quartic_bond");
ff_bond.number_of_parameters = 4;
ff_bond.number_of_parameters = 4;
}
/* ANGLE */

1566
tools/msi2lmp/src/MakeLists.c
View File

File diff suppressed because it is too large Load Diff

342
tools/msi2lmp/src/ReadCarFile.c
View File

@ -1,171 +1,171 @@
/*
* This function opens the .car file and extracts coordinate information
* into the atoms Atom structure
*/
#include "Msi2LMP2.h"
void ReadCarFile(void)
{
char line[MAX_LINE_LENGTH]; /* Stores lines as they are read in */
int k,m,n; /* counters */
int skip; /* lines to skip at beginning of file */
double lowest, highest; /* temp coordinate finding variables */
double total_q;
/* Open .car file for reading */
sprintf(line,"%s.car",rootname);
if (pflag > 0) fprintf(stderr," Reading car file: %s\n",line);
if( (CarF = fopen(line,"r")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
/* Determine Number of molecules & atoms */
rewind(CarF);
no_molecules = -1; /* Set to -1 because counter will be incremented an
extra time at the end of the file */
fgets(line,MAX_LINE_LENGTH,CarF); /* Read header line */
/* Check for periodicity, if present, read cell constants */
if( strncmp(fgets(line,MAX_LINE_LENGTH,CarF),"PBC=ON",6) == 0) {
periodic = 1;
skip = 5; /* Data starts 5 lines from beginning of file */
fgets(line,MAX_LINE_LENGTH,CarF); /* Comment line */
fgets(line,MAX_LINE_LENGTH,CarF); /* Date stamp */
fscanf(CarF,"%*s %lf %lf %lf %lf %lf %lf %*s",
&pbc[0],&pbc[1],&pbc[2],&pbc[3],&pbc[4],&pbc[5]);
if(pbc[3] != 90.0 || pbc[4] != 90.0 || pbc[5] != 90.0) {
fprintf(stderr,"The system is not rectangular- LAMMPS can't handle it!!");
exit(2);
}
}
else {
periodic = 0;
skip = 4;
if (pflag > 1) {
fprintf(stderr," %s is not a periodic system\n", rootname);
fprintf(stderr," Assigning cell parameters based on coordinates\n");
}
fgets(line,MAX_LINE_LENGTH, CarF); /* Comment line */
fgets(line,MAX_LINE_LENGTH, CarF); /* Date Stamp */
}
/* First pass through file -- Count molecules */
while(fgets(line,MAX_LINE_LENGTH,CarF) != NULL )
if( strncmp(line,"end",3) == 0 )
no_molecules++;
/* Allocate space to keep track of the number of atoms within a molecule */
no_atoms = (int *) calloc(no_molecules,sizeof(int));
if ( no_atoms == NULL ) {
fprintf(stderr,"Could not allocate memory for no_atoms\n");
exit(2);
}
/* Second pass through file -- Count atoms */
rewind(CarF);
for(n=0; n < skip; n++) /* Skip beginning lines */
fgets(line,MAX_LINE_LENGTH,CarF);
for(n=0; n < no_molecules; n++)
while( strncmp(fgets(line,MAX_LINE_LENGTH,CarF),"end",3) )
no_atoms[n]++;
for( total_no_atoms=0, n=0; n < no_molecules; n++ )
total_no_atoms += no_atoms[n];
molecule = (struct MoleculeList *) calloc(no_molecules,
sizeof(struct MoleculeList));
if (molecule == NULL) {
fprintf(stderr,"Unable to allocate memory for molecule structure\n");
exit(2);
}
molecule[0].start = 0;
molecule[0].end = no_atoms[0];
for (n=1; n < no_molecules; n++) {
molecule[n].start = molecule[n-1].end;
molecule[n].end = molecule[n].start + no_atoms[n];
}
/* Allocate space for atoms Atom structures */
atoms = (struct Atom *) calloc(total_no_atoms,sizeof(struct Atom));
if( atoms == NULL ) {
fprintf(stderr,"Could not allocate memory for AtomList\n");
exit(2);
}
/* Third pass through file -- Read+Parse Car File */
rewind(CarF);
for(n=0; n < skip; n++)
fgets(line,MAX_LINE_LENGTH,CarF);
for(m=0; m < no_molecules; m++) {
for(k=molecule[m].start; k <
molecule[m].end; k++) {
atoms[k].molecule = m;
atoms[k].no = k;
fscanf(CarF,"%s %lf %lf %lf %*s %*s %s %s %f",
atoms[k].name,
&(atoms[k].x[0]),
&(atoms[k].x[1]),
&(atoms[k].x[2]),
atoms[k].potential,
atoms[k].element,
&(atoms[k].q));
}
fgets(line,MAX_LINE_LENGTH,CarF);
fgets(line,MAX_LINE_LENGTH,CarF);
} /* End m (molecule) loop */
for (total_q=0.0,k=0; k < total_no_atoms; k++)
total_q += atoms[k].q;
if (pflag > 1) {
fprintf(stderr," There are %d atoms in %d molecules in this file\n",
total_no_atoms,no_molecules);
fprintf(stderr," The total charge in the system is %7.3f.\n\n",total_q);
}
/* Search coordinates to find lowest and highest for x, y, and z */
if (periodic == 0) {
for ( k = 0; k < 3; k++) {
lowest = atoms[0].x[k];
highest = atoms[0].x[k];
for ( m = 1; m < total_no_atoms; m++) {
if (atoms[m].x[k] < lowest) lowest = atoms[m].x[k];
if (atoms[m].x[k] > highest) highest = atoms[m].x[k];
}
pbc[k] = lowest;
pbc[k+3] = highest;
}
}
else {
for (k=0; k < 3; k++) {
pbc[k+3] = pbc[k];
pbc[k] = 0.0;
}
}
/* Close .car file */
if (fclose(CarF) !=0) {
fprintf(stderr,"Error closing %s.car\n", rootname);
exit(1);
}
}
/* End ReadCarFile() */
/*
* This function opens the .car file and extracts coordinate information
* into the atoms Atom structure
*/
#include "Msi2LMP2.h"
void ReadCarFile(void)
{
char line[MAX_LINE_LENGTH]; /* Stores lines as they are read in */
int k,m,n; /* counters */
int skip; /* lines to skip at beginning of file */
double lowest, highest; /* temp coordinate finding variables */
double total_q;
/* Open .car file for reading */
sprintf(line,"%s.car",rootname);
if (pflag > 0) fprintf(stderr," Reading car file: %s\n",line);
if( (CarF = fopen(line,"r")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
/* Determine Number of molecules & atoms */
rewind(CarF);
no_molecules = -1; /* Set to -1 because counter will be incremented an
extra time at the end of the file */
fgets(line,MAX_LINE_LENGTH,CarF); /* Read header line */
/* Check for periodicity, if present, read cell constants */
if( strncmp(fgets(line,MAX_LINE_LENGTH,CarF),"PBC=ON",6) == 0) {
periodic = 1;
skip = 5; /* Data starts 5 lines from beginning of file */
fgets(line,MAX_LINE_LENGTH,CarF); /* Comment line */
fgets(line,MAX_LINE_LENGTH,CarF); /* Date stamp */
fscanf(CarF,"%*s %lf %lf %lf %lf %lf %lf %*s",
&pbc[0],&pbc[1],&pbc[2],&pbc[3],&pbc[4],&pbc[5]);
if(pbc[3] != 90.0 || pbc[4] != 90.0 || pbc[5] != 90.0) {
fprintf(stderr,"The system is not rectangular- LAMMPS can't handle it!!");
exit(2);
}
}
else {
periodic = 0;
skip = 4;
if (pflag > 1) {
fprintf(stderr," %s is not a periodic system\n", rootname);
fprintf(stderr," Assigning cell parameters based on coordinates\n");
}
fgets(line,MAX_LINE_LENGTH, CarF); /* Comment line */
fgets(line,MAX_LINE_LENGTH, CarF); /* Date Stamp */
}
/* First pass through file -- Count molecules */
while(fgets(line,MAX_LINE_LENGTH,CarF) != NULL )
if( strncmp(line,"end",3) == 0 )
no_molecules++;
/* Allocate space to keep track of the number of atoms within a molecule */
no_atoms = (int *) calloc(no_molecules,sizeof(int));
if ( no_atoms == NULL ) {
fprintf(stderr,"Could not allocate memory for no_atoms\n");
exit(2);
}
/* Second pass through file -- Count atoms */
rewind(CarF);
for(n=0; n < skip; n++) /* Skip beginning lines */
fgets(line,MAX_LINE_LENGTH,CarF);
for(n=0; n < no_molecules; n++)
while( strncmp(fgets(line,MAX_LINE_LENGTH,CarF),"end",3) )
no_atoms[n]++;
for( total_no_atoms=0, n=0; n < no_molecules; n++ )
total_no_atoms += no_atoms[n];
molecule = (struct MoleculeList *) calloc(no_molecules,
sizeof(struct MoleculeList));
if (molecule == NULL) {
fprintf(stderr,"Unable to allocate memory for molecule structure\n");
exit(2);
}
molecule[0].start = 0;
molecule[0].end = no_atoms[0];
for (n=1; n < no_molecules; n++) {
molecule[n].start = molecule[n-1].end;
molecule[n].end = molecule[n].start + no_atoms[n];
}
/* Allocate space for atoms Atom structures */
atoms = (struct Atom *) calloc(total_no_atoms,sizeof(struct Atom));
if( atoms == NULL ) {
fprintf(stderr,"Could not allocate memory for AtomList\n");
exit(2);
}
/* Third pass through file -- Read+Parse Car File */
rewind(CarF);
for(n=0; n < skip; n++)
fgets(line,MAX_LINE_LENGTH,CarF);
for(m=0; m < no_molecules; m++) {
for(k=molecule[m].start; k <
molecule[m].end; k++) {
atoms[k].molecule = m;
atoms[k].no = k;
fscanf(CarF,"%s %lf %lf %lf %*s %*s %s %s %f",
atoms[k].name,
&(atoms[k].x[0]),
&(atoms[k].x[1]),
&(atoms[k].x[2]),
atoms[k].potential,
atoms[k].element,
&(atoms[k].q));
}
fgets(line,MAX_LINE_LENGTH,CarF);
fgets(line,MAX_LINE_LENGTH,CarF);
} /* End m (molecule) loop */
for (total_q=0.0,k=0; k < total_no_atoms; k++)
total_q += atoms[k].q;
if (pflag > 1) {
fprintf(stderr," There are %d atoms in %d molecules in this file\n",
total_no_atoms,no_molecules);
fprintf(stderr," The total charge in the system is %7.3f.\n\n",total_q);
}
/* Search coordinates to find lowest and highest for x, y, and z */
if (periodic == 0) {
for ( k = 0; k < 3; k++) {
lowest = atoms[0].x[k];
highest = atoms[0].x[k];
for ( m = 1; m < total_no_atoms; m++) {
if (atoms[m].x[k] < lowest) lowest = atoms[m].x[k];
if (atoms[m].x[k] > highest) highest = atoms[m].x[k];
}
pbc[k] = lowest;
pbc[k+3] = highest;
}
}
else {
for (k=0; k < 3; k++) {
pbc[k+3] = pbc[k];
pbc[k] = 0.0;
}
}
/* Close .car file */
if (fclose(CarF) !=0) {
fprintf(stderr,"Error closing %s.car\n", rootname);
exit(1);
}
}
/* End ReadCarFile() */

850
tools/msi2lmp/src/ReadMdfFile.c
View File

@ -1,425 +1,425 @@
/******************************
*
* This function opens the .mdf file and extracts connectivity information
* into the atoms Atom structure. It also updates the charge from the .car
* file because the charge in the .mdf file has more significant figures.
*
*/
#include "Msi2LMP2.h"
/* Prototype for function to process a single atom
Returns int that flags end of data file */
int get_molecule(char line[], int connect_col_no,
int q_col_no, int *counter);
/* Prototype for function that takes connectivty record as stated in
.mdf file and fills in any default values */
void MakeConnectFullForm(int *counter);
/* prototype for function to clean strange characters out of strings */
void clean_string(char *);
void ReadMdfFile(void)
{
char line[MAX_LINE_LENGTH]; /* Temporary storage for reading lines */
char *col_no; /* Pointer to column number stored as char */
char *col_name; /* Pointer to column name */
int connect_col_no = 0; /* Column number where connection info begins */
int q_col_no = 0; /* Column number containg charge information */
int column_flag=0; /* Flag for finding connect and q col no */
int atom_counter=0; /* Keeps track of current atom number */
int p_flag = 1; /* return value from ProcessConnections() */
int i,j,k,kk,l,n,match,match2,status;
char *temp_string;
char *temp_residue;
char *temp_atom_name;
char *sptr;
char *molecule_name;
unsigned char at_end = 0;
/* Open .mdf file for reading */
sprintf(line,"%s.mdf",rootname);
if (pflag > 0) fprintf(stderr," Reading mdf file: %s\n",line);
if ((MdfF = fopen(line,"r")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
while (!at_end) {
sptr = fgets(line,MAX_LINE_LENGTH,MdfF);
if (sptr != NULL) {
clean_string(line);
if (strncmp(line,"#end",4) == 0) {
at_end = 1;
}
else if (strncmp(line,"@column",7) == 0) {
temp_string = strtok(line," ");
col_no = strtok(NULL," ");
col_name = strtok(NULL," ");
if (strncmp(col_name,"charge",6) == 0) {
if (strlen(col_name) < 8) {
q_col_no = atoi(col_no);
}
}
else if (strncmp(col_name,"connect",7) == 0) {
connect_col_no = atoi(col_no);
}
}
else if (strncmp(line,"@molecule",9) == 0) {
temp_string = strtok(line," ");
molecule_name = strtok(NULL," ");
if ((q_col_no == 0) | (connect_col_no == 0)) {
fprintf(stderr,"Unable to process molecule without knowing charge\n");
fprintf(stderr,"and connections columns\n");
}
sptr = fgets(line,MAX_LINE_LENGTH,MdfF);
status = get_molecule(line,connect_col_no,q_col_no,&atom_counter);
if (status == 0) {
fprintf(stderr,"Trouble reading molecule - exiting\n");
}
}
else {
}
}
else {
fprintf(stderr,"End of File found or error reading line\n");
at_end = 1;
}
}
/* Next build list of residues for each molecule This will
facilitate assigning connections numbers as well as figuring
out bonds, angles, etc. This first loop just figures out the
number of residues in each molecule and allocates memory to
store information for each residue. The second loop fills
in starting and ending atom positions for each residue
*/
temp_string = calloc(16,sizeof(char));
temp_string[15] = '\0';
for (n=0; n < no_molecules; n++) {
molecule[n].no_residues = 1;
strncpy(temp_string,atoms[molecule[n].start].residue_string,16);
for (i=molecule[n].start+1; i < molecule[n].end; i++) {
if (strncmp(temp_string,atoms[i].residue_string,16) != 0) {
molecule[n].no_residues++;
strncpy(temp_string,atoms[i].residue_string,16);
}
}
molecule[n].residue = (struct ResidueList *)
calloc(molecule[n].no_residues, sizeof(struct ResidueList));
if (molecule[n].residue == NULL) {
fprintf(stderr,"Unable to allocate memory for residue list - molecule %d\n",
n);
exit(1);
}
}
for (n=0; n < no_molecules; n++) {
j = 0;
strncpy(molecule[n].residue[j].name,
atoms[molecule[n].start].residue_string,16);
molecule[n].residue[j].start = molecule[n].start;
for (i=molecule[n].start+1; i < molecule[n].end; i++) {
if (strncmp(molecule[n].residue[j].name,
atoms[i].residue_string,16) != 0) {
molecule[n].residue[j].end = i;
molecule[n].residue[++j].start = i;
strncpy(molecule[n].residue[j].name,atoms[i].residue_string,16);
}
}
molecule[n].residue[j].end = molecule[n].end;
/*
fprintf(stderr,"Molecule %d has %d residues",n,molecule[n].no_residues);
for (i=0; i < molecule[n].no_residues; i++) {
fprintf(stderr," %s",molecule[n].residue[i].name);
}
fprintf(stderr,"\n");
for (i=molecule[n].start; i < molecule[n].end; i++) {
fprintf(stderr," atom %d residue %s\n",i,atoms[i].residue_string);
}
fprintf(stderr," residue %s start %d end %d\n",molecule[n].residue[i].name,
molecule[n].residue[i].start,molecule[n].residue[i].end);
}
*/
}
/* Assign atom names in connections[] to corresponding atom numbers */
for (n=0; n < no_molecules; n++) {
for (j=0; j < molecule[n].no_residues; j++) {
for (i=molecule[n].residue[j].start; i < molecule[n].residue[j].end;
i++) {
for (l=0; l < atoms[i].no_connect; l++) {
strncpy(temp_string,atoms[i].connections[l],16);
temp_residue = strtok(temp_string,":");
temp_atom_name = strtok(NULL,"%");
if (strcmp(temp_residue,molecule[n].residue[j].name) == 0) {
/* atom and connection are part of same residue
Search names on just that residue */
k = molecule[n].residue[j].start;
match = 0;
while (!match && (k < molecule[n].residue[j].end)) {
if (strcmp(atoms[k].name,temp_atom_name) == 0) {
atoms[i].conn_no[l] = k;
match = 1;
}
else
k++;
}
if (match == 0) {
fprintf(stderr,"Unable to resolve atom number of atom %d conn %d string %s:%s\n Something is wrong in the MDF file\n",
i,l,temp_residue,temp_atom_name);
exit(1);
}
}
else {
/* atom and connection are on different residues
First find the residue that the connection is
on then loop over its atoms
*/
k=0;
match = 0;
while (!match && (k < molecule[n].no_residues)) {
if (strcmp(temp_residue,molecule[n].residue[k].name) == 0) {
kk = molecule[n].residue[k].start;
match2 = 0;
while (!match2 && (kk < molecule[n].residue[k].end)) {
if (strcmp(atoms[kk].name,temp_atom_name) == 0) {
atoms[i].conn_no[l] = kk;
match2 = 1;
}
else
kk++;
}
if (match2 == 0) {
fprintf(stderr,"Unable to resolve atom number of atom %d conn %d string %s\n Something is wrong in the MDF file\n",
i,l,atoms[i].connections[l]);
exit(1);
}
match = 1;
}
else
k++;
}
if (match == 0) {
fprintf(stderr,"Unable to find residue associated with conn %d %s on atom %d\n Something is wrong in the MDF file\n", l,atoms[i].connections[l],i);
exit(1);
}
} /* end if */
} /* l - loop over connections on atom i */
} /* i - loop on atoms in residue j molecule n */
} /* j - loop on residues in molecule n */
} /* n - loop over molecules */
free(temp_string);
/*
for (n=0; n < no_molecules; n++) {
fprintf(stderr,"Molecule %d has %d residues\n",n,molecule[n].no_residues);
for (j=0; j < molecule[n].no_residues; j++) {
fprintf(stderr," Residue %d named %s\n",j,molecule[n].residue[j].name);
for (i=molecule[n].residue[j].start; i < molecule[n].residue[j].end;
i++) {
fprintf(stderr," Atom %d type %s connected to ",i,atoms[i].potential);
for (l=0; l < atoms[i].no_connect; l++) fprintf(stderr," %d ",
atoms[i].conn_no[l]);
fprintf(stderr,"\n");
}
}
}
*/
/* Close .mdf file */
if (fclose(MdfF) !=0) {
printf("Error closing %s.car\n", rootname);
exit(1);
}
} /* End ReadMdfFile function */
/*--------------------- get_molecule Function-----------------------*/
int get_molecule(char *line, int connect_col_no, int q_col_no,
int *counter)
{
char *cur_field; /* For storing current string token */
int i; /* Used in loop counters */
int connect_no; /* Connection number within atom */
int r_val = 1; /* Return value. 1 = successful
0 = EOF encountered */
/* Loop over atoms */
/* blank line signals end of molecule*/
while(!blank_line(fgets(line,MAX_LINE_LENGTH,MdfF))) {
/* while(strlen(fgets(line,MAX_LINE_LENGTH,MdfF)) > 2) { */
clean_string(line);
/* Get atom name */
cur_field = strtok(line,":");
sscanf(cur_field, "%s", atoms[*counter].residue_string);
cur_field = strtok(NULL," ");
/* Compare atom name with that in .car file */
if (strcmp(atoms[*counter].name, cur_field)) {
fprintf(stderr,"Names %s from .car file and %s from .mdf file do not match\n",
atoms[*counter].name, cur_field);
fprintf(stderr,"counter = %d\n",*counter);
fprintf(stderr,"Program Terminating\n");
exit(4);
}
/* Skip unwanted fields until charge column, then update charge */
for (i=1; i < q_col_no; i++) strtok(NULL," ");
cur_field = strtok(NULL, " ");
atoms[*counter].q = atof(cur_field);
/* Continue skipping unwanted fields until connectivity records begin */
for ( i = (q_col_no + 1); i < connect_col_no; i++) strtok(NULL," ");
/* Process connections */
connect_no = 0; /* reset connections counter */
while ((cur_field = strtok(NULL," ")) && (connect_no < MAX_CONNECTIONS)) {
sscanf(cur_field, "%s", atoms[*counter].connections[connect_no++]);
}
atoms[*counter].no_connect = connect_no;
MakeConnectFullForm(counter);
(*counter)++;
} /* End atom processing loop */
return r_val;
} /* End get_molecule function */
/*------------------------MakeConnectFullForm Function--------------------*/
void MakeConnectFullForm(int *counter) {
/* This function processes the connection names after all connections
for an atom have been read in.
It replaces any short forms that use implied default values
with the full form connectivity record */
int i; /* Counter for character array */
int j; /* loop counter */
char tempname[MAX_STRING]; /* name of connection */
char tempcell[10]; /* Values from connectivity record */
char tempsym[5]; /* " " */
char tempbo[6]; /* " " */
char *charptr;
for ( j = 0; j < atoms[*counter].no_connect; j++) {
/* If not full name, make name full */
if (strchr(atoms[*counter].connections[j],':') == NULL) {
strcpy(tempname,atoms[*counter].residue_string);
strcat(tempname,":");
strcat(tempname,
atoms[*counter].connections[j]);
sscanf(tempname, "%s",
atoms[*counter].connections[j]);
}
else
sscanf(atoms[*counter].connections[j], "%s", tempname);
/* Set cell variables */
i=0;
charptr = (strchr(tempname,'%'));
if (charptr != NULL) {
while ( *charptr!='#' && *charptr!='/' && *charptr!='\000')
tempcell[i++] = *(charptr++);
tempcell[i] = '\000';
}
else
strcpy(tempcell, "%000");
/* Set symmetry variables
-- If not 1, cannot handle at this time */
i = 0;
charptr = (strchr(tempname,'#'));
if (charptr != NULL) {
while (*charptr != '/' && *charptr !='\000') {
tempsym[i++] = *(charptr++);
if ((i==2) && (tempsym[1] != '1')) {
fprintf(stderr,"Msi2LMP is not equipped to handle symmetry operations\n");
exit(5);
}
}
tempsym[i] = '\000';
}
else
strcpy(tempsym, "#1");
/* Set bond order and record in data structure */
i = 0;
charptr = strchr(tempname,'/');
if (charptr != NULL) {
charptr++;
while (*charptr != '\000')
tempbo[i++] = *(charptr++);
tempbo[i] = '\000';
}
else
strcpy(tempbo, "1.0");
atoms[*counter].bond_order[j] = atof(tempbo);
/* Build connection name and store in atoms data structure */
strtok( tempname, "%#/");
strcat( tempname, tempcell);
strcat( tempname, tempsym);
strcat( tempname, "/");
strcat( tempname, tempbo);
if (strlen(tempname) > 25) fprintf(stderr,"tempname overrun %s\n",tempname);
sscanf( tempname, "%s", atoms[*counter].connections[j]);
}/*End for loop*/
}/* End function MakeNameLong
*/
void clean_string(char *string) {
int i,n;
short k;
n = strlen(string);
for (i=0; i < n; i++) {
k = (short)string[i];
if ((k<32) | (k>127)) string[i] = '\0';
}
}
/******************************
*
* This function opens the .mdf file and extracts connectivity information
* into the atoms Atom structure. It also updates the charge from the .car
* file because the charge in the .mdf file has more significant figures.
*
*/
#include "Msi2LMP2.h"
/* Prototype for function to process a single atom
Returns int that flags end of data file */
int get_molecule(char line[], int connect_col_no,
int q_col_no, int *counter);
/* Prototype for function that takes connectivty record as stated in
.mdf file and fills in any default values */
void MakeConnectFullForm(int *counter);
/* prototype for function to clean strange characters out of strings */
void clean_string(char *);
void ReadMdfFile(void)
{
char line[MAX_LINE_LENGTH]; /* Temporary storage for reading lines */
char *col_no; /* Pointer to column number stored as char */
char *col_name; /* Pointer to column name */
int connect_col_no = 0; /* Column number where connection info begins */
int q_col_no = 0; /* Column number containg charge information */
int column_flag=0; /* Flag for finding connect and q col no */
int atom_counter=0; /* Keeps track of current atom number */
int p_flag = 1; /* return value from ProcessConnections() */
int i,j,k,kk,l,n,match,match2,status;
char *temp_string;
char *temp_residue;
char *temp_atom_name;
char *sptr;
char *molecule_name;
unsigned char at_end = 0;
/* Open .mdf file for reading */
sprintf(line,"%s.mdf",rootname);
if (pflag > 0) fprintf(stderr," Reading mdf file: %s\n",line);
if ((MdfF = fopen(line,"r")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
while (!at_end) {
sptr = fgets(line,MAX_LINE_LENGTH,MdfF);
if (sptr != NULL) {
clean_string(line);
if (strncmp(line,"#end",4) == 0) {
at_end = 1;
}
else if (strncmp(line,"@column",7) == 0) {
temp_string = strtok(line," ");
col_no = strtok(NULL," ");
col_name = strtok(NULL," ");
if (strncmp(col_name,"charge",6) == 0) {
if (strlen(col_name) < 8) {
q_col_no = atoi(col_no);
}
}
else if (strncmp(col_name,"connect",7) == 0) {
connect_col_no = atoi(col_no);
}
}
else if (strncmp(line,"@molecule",9) == 0) {
temp_string = strtok(line," ");
molecule_name = strtok(NULL," ");
if ((q_col_no == 0) | (connect_col_no == 0)) {
fprintf(stderr,"Unable to process molecule without knowing charge\n");
fprintf(stderr,"and connections columns\n");
}
sptr = fgets(line,MAX_LINE_LENGTH,MdfF);
status = get_molecule(line,connect_col_no,q_col_no,&atom_counter);
if (status == 0) {
fprintf(stderr,"Trouble reading molecule - exiting\n");
}
}
else {
}
}
else {
fprintf(stderr,"End of File found or error reading line\n");
at_end = 1;
}
}
/* Next build list of residues for each molecule This will
facilitate assigning connections numbers as well as figuring
out bonds, angles, etc. This first loop just figures out the
number of residues in each molecule and allocates memory to
store information for each residue. The second loop fills
in starting and ending atom positions for each residue
*/
temp_string = calloc(16,sizeof(char));
temp_string[15] = '\0';
for (n=0; n < no_molecules; n++) {
molecule[n].no_residues = 1;
strncpy(temp_string,atoms[molecule[n].start].residue_string,16);
for (i=molecule[n].start+1; i < molecule[n].end; i++) {
if (strncmp(temp_string,atoms[i].residue_string,16) != 0) {
molecule[n].no_residues++;
strncpy(temp_string,atoms[i].residue_string,16);
}
}
molecule[n].residue = (struct ResidueList *)
calloc(molecule[n].no_residues, sizeof(struct ResidueList));
if (molecule[n].residue == NULL) {
fprintf(stderr,"Unable to allocate memory for residue list - molecule %d\n",
n);
exit(1);
}
}
for (n=0; n < no_molecules; n++) {
j = 0;
strncpy(molecule[n].residue[j].name,
atoms[molecule[n].start].residue_string,16);
molecule[n].residue[j].start = molecule[n].start;
for (i=molecule[n].start+1; i < molecule[n].end; i++) {
if (strncmp(molecule[n].residue[j].name,
atoms[i].residue_string,16) != 0) {
molecule[n].residue[j].end = i;
molecule[n].residue[++j].start = i;
strncpy(molecule[n].residue[j].name,atoms[i].residue_string,16);
}
}
molecule[n].residue[j].end = molecule[n].end;
/*
fprintf(stderr,"Molecule %d has %d residues",n,molecule[n].no_residues);
for (i=0; i < molecule[n].no_residues; i++) {
fprintf(stderr," %s",molecule[n].residue[i].name);
}
fprintf(stderr,"\n");
for (i=molecule[n].start; i < molecule[n].end; i++) {
fprintf(stderr," atom %d residue %s\n",i,atoms[i].residue_string);
}
fprintf(stderr," residue %s start %d end %d\n",molecule[n].residue[i].name,
molecule[n].residue[i].start,molecule[n].residue[i].end);
}
*/
}
/* Assign atom names in connections[] to corresponding atom numbers */
for (n=0; n < no_molecules; n++) {
for (j=0; j < molecule[n].no_residues; j++) {
for (i=molecule[n].residue[j].start; i < molecule[n].residue[j].end;
i++) {
for (l=0; l < atoms[i].no_connect; l++) {
strncpy(temp_string,atoms[i].connections[l],16);
temp_residue = strtok(temp_string,":");
temp_atom_name = strtok(NULL,"%");
if (strcmp(temp_residue,molecule[n].residue[j].name) == 0) {
/* atom and connection are part of same residue
Search names on just that residue */
k = molecule[n].residue[j].start;
match = 0;
while (!match && (k < molecule[n].residue[j].end)) {
if (strcmp(atoms[k].name,temp_atom_name) == 0) {
atoms[i].conn_no[l] = k;
match = 1;
}
else
k++;
}
if (match == 0) {
fprintf(stderr,"Unable to resolve atom number of atom %d conn %d string %s:%s\n Something is wrong in the MDF file\n",
i,l,temp_residue,temp_atom_name);
exit(1);
}
}
else {
/* atom and connection are on different residues
First find the residue that the connection is
on then loop over its atoms
*/
k=0;
match = 0;
while (!match && (k < molecule[n].no_residues)) {
if (strcmp(temp_residue,molecule[n].residue[k].name) == 0) {
kk = molecule[n].residue[k].start;
match2 = 0;
while (!match2 && (kk < molecule[n].residue[k].end)) {
if (strcmp(atoms[kk].name,temp_atom_name) == 0) {
atoms[i].conn_no[l] = kk;
match2 = 1;
}
else
kk++;
}
if (match2 == 0) {
fprintf(stderr,"Unable to resolve atom number of atom %d conn %d string %s\n Something is wrong in the MDF file\n",
i,l,atoms[i].connections[l]);
exit(1);
}
match = 1;
}
else
k++;
}
if (match == 0) {
fprintf(stderr,"Unable to find residue associated with conn %d %s on atom %d\n Something is wrong in the MDF file\n", l,atoms[i].connections[l],i);
exit(1);
}
} /* end if */
} /* l - loop over connections on atom i */
} /* i - loop on atoms in residue j molecule n */
} /* j - loop on residues in molecule n */
} /* n - loop over molecules */
free(temp_string);
/*
for (n=0; n < no_molecules; n++) {
fprintf(stderr,"Molecule %d has %d residues\n",n,molecule[n].no_residues);
for (j=0; j < molecule[n].no_residues; j++) {
fprintf(stderr," Residue %d named %s\n",j,molecule[n].residue[j].name);
for (i=molecule[n].residue[j].start; i < molecule[n].residue[j].end;
i++) {
fprintf(stderr," Atom %d type %s connected to ",i,atoms[i].potential);
for (l=0; l < atoms[i].no_connect; l++) fprintf(stderr," %d ",
atoms[i].conn_no[l]);
fprintf(stderr,"\n");
}
}
}
*/
/* Close .mdf file */
if (fclose(MdfF) !=0) {
printf("Error closing %s.car\n", rootname);
exit(1);
}
} /* End ReadMdfFile function */
/*--------------------- get_molecule Function-----------------------*/
int get_molecule(char *line, int connect_col_no, int q_col_no,
int *counter)
{
char *cur_field; /* For storing current string token */
int i; /* Used in loop counters */
int connect_no; /* Connection number within atom */
int r_val = 1; /* Return value. 1 = successful
0 = EOF encountered */
/* Loop over atoms */
/* blank line signals end of molecule*/
while(!blank_line(fgets(line,MAX_LINE_LENGTH,MdfF))) {
/* while(strlen(fgets(line,MAX_LINE_LENGTH,MdfF)) > 2) { */
clean_string(line);
/* Get atom name */
cur_field = strtok(line,":");
sscanf(cur_field, "%s", atoms[*counter].residue_string);
cur_field = strtok(NULL," ");
/* Compare atom name with that in .car file */
if (strcmp(atoms[*counter].name, cur_field)) {
fprintf(stderr,"Names %s from .car file and %s from .mdf file do not match\n",
atoms[*counter].name, cur_field);
fprintf(stderr,"counter = %d\n",*counter);
fprintf(stderr,"Program Terminating\n");
exit(4);
}
/* Skip unwanted fields until charge column, then update charge */
for (i=1; i < q_col_no; i++) strtok(NULL," ");
cur_field = strtok(NULL, " ");
atoms[*counter].q = atof(cur_field);
/* Continue skipping unwanted fields until connectivity records begin */
for ( i = (q_col_no + 1); i < connect_col_no; i++) strtok(NULL," ");
/* Process connections */
connect_no = 0; /* reset connections counter */
while ((cur_field = strtok(NULL," ")) && (connect_no < MAX_CONNECTIONS)) {
sscanf(cur_field, "%s", atoms[*counter].connections[connect_no++]);
}
atoms[*counter].no_connect = connect_no;
MakeConnectFullForm(counter);
(*counter)++;
} /* End atom processing loop */
return r_val;
} /* End get_molecule function */
/*------------------------MakeConnectFullForm Function--------------------*/
void MakeConnectFullForm(int *counter) {
/* This function processes the connection names after all connections
for an atom have been read in.
It replaces any short forms that use implied default values
with the full form connectivity record */
int i; /* Counter for character array */
int j; /* loop counter */
char tempname[MAX_STRING]; /* name of connection */
char tempcell[10]; /* Values from connectivity record */
char tempsym[5]; /* " " */
char tempbo[6]; /* " " */
char *charptr;
for ( j = 0; j < atoms[*counter].no_connect; j++) {
/* If not full name, make name full */
if (strchr(atoms[*counter].connections[j],':') == NULL) {
strcpy(tempname,atoms[*counter].residue_string);
strcat(tempname,":");
strcat(tempname,
atoms[*counter].connections[j]);
sscanf(tempname, "%s",
atoms[*counter].connections[j]);
}
else
sscanf(atoms[*counter].connections[j], "%s", tempname);
/* Set cell variables */
i=0;
charptr = (strchr(tempname,'%'));
if (charptr != NULL) {
while ( *charptr!='#' && *charptr!='/' && *charptr!='\000')
tempcell[i++] = *(charptr++);
tempcell[i] = '\000';
}
else
strcpy(tempcell, "%000");
/* Set symmetry variables
-- If not 1, cannot handle at this time */
i = 0;
charptr = (strchr(tempname,'#'));
if (charptr != NULL) {
while (*charptr != '/' && *charptr !='\000') {
tempsym[i++] = *(charptr++);
if ((i==2) && (tempsym[1] != '1')) {
fprintf(stderr,"Msi2LMP is not equipped to handle symmetry operations\n");
exit(5);
}
}
tempsym[i] = '\000';
}
else
strcpy(tempsym, "#1");
/* Set bond order and record in data structure */
i = 0;
charptr = strchr(tempname,'/');
if (charptr != NULL) {
charptr++;
while (*charptr != '\000')
tempbo[i++] = *(charptr++);
tempbo[i] = '\000';
}
else
strcpy(tempbo, "1.0");
atoms[*counter].bond_order[j] = atof(tempbo);
/* Build connection name and store in atoms data structure */
strtok( tempname, "%#/");
strcat( tempname, tempcell);
strcat( tempname, tempsym);
strcat( tempname, "/");
strcat( tempname, tempbo);
if (strlen(tempname) > 25) fprintf(stderr,"tempname overrun %s\n",tempname);
sscanf( tempname, "%s", atoms[*counter].connections[j]);
}/*End for loop*/
}/* End function MakeNameLong
*/
void clean_string(char *string) {
int i,n;
short k;
n = strlen(string);
for (i=0; i < n; i++) {
k = (short)string[i];
if ((k<32) | (k>127)) string[i] = '\0';
}
}

436
tools/msi2lmp/src/SearchAndFill.c
View File

@ -1,218 +1,218 @@
/****************************
*
* This function first allocates memory to the forcefield item
* structures and then reads parameters from the forcefield file into the
* allocated memory
*
*/
#include "Forcefield.h"
#include "Msi2LMP2.h"
unsigned char string_match(char *,char *);
void SearchAndFill(struct FrcFieldItem *item)
{
int i,j; /* counters */
int got_it = 0;
int ctr = 0;
long file_pos;
char line[MAX_LINE] = "empty";
char *charptr,*status;
extern FILE *FrcF;
/***********************ALLOCATE MEMORY FOR STRUCTURE ********************/
/* Read and discard lines until keyword is found */
rewind(FrcF);
while ((got_it == 0)) {
status = fgets( line, MAX_LINE, FrcF );
if (status == NULL) {
fprintf(stderr," Unable to find forcefield keyword %s\n",item->keyword);
fprintf(stderr," Check consistency of forcefield name and class \n");
fprintf(stderr," Exiting....\n");
exit(1);
}
if (line[0] == '#') {
if (string_match(strtok(line," '\t'("),item->keyword)) got_it = 1;
}
/* if (strncmp(line, item->keyword,strlen(item->keyword))==0) got_it = 1; */
}
file_pos = ftell(FrcF);
/* Count the number of lines until next item is found */
while( strncmp(fgets(line,MAX_LINE,FrcF), "#", 1) != 0 )
ctr++;
/* Allocate the memory using calloc */
item->data = calloc(ctr, sizeof(struct FrcFieldData));
if (item->data == NULL) {
fprintf(stderr,"Could not allocate memory to %s\n", item->keyword);
exit(2);
}
/********************FILL PARAMETERS AND EQUIVALENCES ********************/
/* Read lines until keyword is found */
fseek(FrcF,file_pos,SEEK_SET);
strcpy(line,"empty");
/* Read lines until data starts (when !--- is found) */
ctr = 0;
while ( strncmp(line,"!---", 4) != 0 ) {
fgets(line, MAX_LINE, FrcF);
}
/* Get first line of data that isn't commented out */
fgets(line, MAX_LINE, FrcF);
while (strncmp(line,"!",1) == 0) {
fgets( line, MAX_LINE, FrcF);
}
/* Read data into structure */
while( strncmp( line, "#", 1 ) != 0 ) {
float version;
int reference,replace;
char atom_types[5][5];
double parameters[8];
/* version number and reference number */
version = atof(strtok(line, " "));
reference = atoi(strtok(NULL, " "));
/* equivalences */
for(i = 0; i < item->number_of_members; i++ ) {
sscanf(strtok(NULL, " "), "%s", atom_types[i]);
}
/* parameters -- Because of symmetrical terms, bonang, angtor, and
endbontor have to be treated carefully */
for( i = 0; i < item->number_of_parameters; i++ ) {
charptr = strtok(NULL, " ");
if(charptr == NULL) {
for ( j = i; j < item->number_of_parameters; j++ )
parameters[j] = parameters[j-i];
break;
}
else {
parameters[i] = atof(charptr);
}
}
/* Search for matching sets of atom types.
If found and the version number is greater, substitute
the current set of parameters in place of the found set.
Otherwise, add the current set of parameters to the
list.
*/
replace = ctr;
for (j=0; j < ctr; j++) {
int k=0;
int match = 1;
while (match && (k < item->number_of_members)) {
if (strncmp(item->data[j].ff_types[k],atom_types[k],5) == 0)
k++;
else
match = 0;
}
if (match == 1) {
replace = j;
break;
}
}
if (replace != ctr) {
if (version > item->data[replace].ver) {
if (pflag > 1) {
fprintf(stderr," Using higher version of parameters for");
fprintf(stderr," %s ",item->keyword);
for (i=0; i < item->number_of_members; i++)
fprintf(stderr,"%s ",atom_types[i]);
fprintf(stderr," version %3.2f\n",version);
}
item->data[replace].ver = version;
item->data[replace].ref = reference;
for (i=0; i < item->number_of_members; i++) {
strncpy(item->data[replace].ff_types[i],atom_types[i],5);
}
for (i=0; i < item->number_of_parameters; i++) {
item->data[replace].ff_param[i] = parameters[i];
}
}
else {
if (pflag > 1) {
fprintf(stderr," Using higher version of parameters for");
fprintf(stderr," %s ",item->keyword);
for (i=0; i < item->number_of_members; i++)
fprintf(stderr,"%s ",item->data[replace].ff_types[i]);
fprintf(stderr," version %3.2f\n",item->data[replace].ver);
}
}
}
else {
item->data[ctr].ver = version;
item->data[ctr].ref = reference;
for (i=0; i < item->number_of_members; i++) {
strncpy(item->data[ctr].ff_types[i],atom_types[i],5);
}
for (i=0; i < item->number_of_parameters; i++) {
item->data[ctr].ff_param[i] = parameters[i];
}
ctr++;
}
fgets( line, MAX_LINE, FrcF);
/*if blank line encountered, get next */
while((blank_line(line)) ||
(strncmp(line,"!",1) == 0)) {
fgets( line, MAX_LINE, FrcF);
}
}
item->entries = ctr;
/*Debugging
fprintf(stderr,"\n%s\n", item->keyword);
for(i=0;i<ctr;i++) {
for(j=0;j<item->number_of_members;j++)
fprintf(stderr,"%3s ", item->data[i].ff_equiv[j]);
fprintf(stderr," ");
for(j=0;j<item->number_of_parameters;j++)
fprintf(stderr,"%10.5f ",item->data[i].ff_param[j]);
fprintf(stderr,"\n");
}
*/
}
unsigned char string_match(char *string1,char *string2)
{
int len1,len2;
len1 = strlen(string1);
len2 = strlen(string2);
if (len1 != len2) {
return 0;
}
else {
if (strncmp(string1,string2,len1) == 0) {
return 1;
}
else {
return 0;
}
}
}
/****************************
*
* This function first allocates memory to the forcefield item
* structures and then reads parameters from the forcefield file into the
* allocated memory
*
*/
#include "Forcefield.h"
#include "Msi2LMP2.h"
unsigned char string_match(char *,char *);
void SearchAndFill(struct FrcFieldItem *item)
{
int i,j; /* counters */
int got_it = 0;
int ctr = 0;
long file_pos;
char line[MAX_LINE] = "empty";
char *charptr,*status;
extern FILE *FrcF;
/***********************ALLOCATE MEMORY FOR STRUCTURE ********************/
/* Read and discard lines until keyword is found */
rewind(FrcF);
while ((got_it == 0)) {
status = fgets( line, MAX_LINE, FrcF );
if (status == NULL) {
fprintf(stderr," Unable to find forcefield keyword %s\n",item->keyword);
fprintf(stderr," Check consistency of forcefield name and class \n");
fprintf(stderr," Exiting....\n");
exit(1);
}
if (line[0] == '#') {
if (string_match(strtok(line," '\t'("),item->keyword)) got_it = 1;
}
/* if (strncmp(line, item->keyword,strlen(item->keyword))==0) got_it = 1; */
}
file_pos = ftell(FrcF);
/* Count the number of lines until next item is found */
while( strncmp(fgets(line,MAX_LINE,FrcF), "#", 1) != 0 )
ctr++;
/* Allocate the memory using calloc */
item->data = calloc(ctr, sizeof(struct FrcFieldData));
if (item->data == NULL) {
fprintf(stderr,"Could not allocate memory to %s\n", item->keyword);
exit(2);
}
/********************FILL PARAMETERS AND EQUIVALENCES ********************/
/* Read lines until keyword is found */
fseek(FrcF,file_pos,SEEK_SET);
strcpy(line,"empty");
/* Read lines until data starts (when !--- is found) */
ctr = 0;
while ( strncmp(line,"!---", 4) != 0 ) {
fgets(line, MAX_LINE, FrcF);
}
/* Get first line of data that isn't commented out */
fgets(line, MAX_LINE, FrcF);
while (strncmp(line,"!",1) == 0) {
fgets( line, MAX_LINE, FrcF);
}
/* Read data into structure */
while( strncmp( line, "#", 1 ) != 0 ) {
float version;
int reference,replace;
char atom_types[5][5];
double parameters[8];
/* version number and reference number */
version = atof(strtok(line, " "));
reference = atoi(strtok(NULL, " "));
/* equivalences */
for(i = 0; i < item->number_of_members; i++ ) {
sscanf(strtok(NULL, " "), "%s", atom_types[i]);
}
/* parameters -- Because of symmetrical terms, bonang, angtor, and
endbontor have to be treated carefully */
for( i = 0; i < item->number_of_parameters; i++ ) {
charptr = strtok(NULL, " ");
if(charptr == NULL) {
for ( j = i; j < item->number_of_parameters; j++ )
parameters[j] = parameters[j-i];
break;
}
else {
parameters[i] = atof(charptr);
}
}
/* Search for matching sets of atom types.
If found and the version number is greater, substitute
the current set of parameters in place of the found set.
Otherwise, add the current set of parameters to the
list.
*/
replace = ctr;
for (j=0; j < ctr; j++) {
int k=0;
int match = 1;
while (match && (k < item->number_of_members)) {
if (strncmp(item->data[j].ff_types[k],atom_types[k],5) == 0)
k++;
else
match = 0;
}
if (match == 1) {
replace = j;
break;
}
}
if (replace != ctr) {
if (version > item->data[replace].ver) {
if (pflag > 1) {
fprintf(stderr," Using higher version of parameters for");
fprintf(stderr," %s ",item->keyword);
for (i=0; i < item->number_of_members; i++)
fprintf(stderr,"%s ",atom_types[i]);
fprintf(stderr," version %3.2f\n",version);
}
item->data[replace].ver = version;
item->data[replace].ref = reference;
for (i=0; i < item->number_of_members; i++) {
strncpy(item->data[replace].ff_types[i],atom_types[i],5);
}
for (i=0; i < item->number_of_parameters; i++) {
item->data[replace].ff_param[i] = parameters[i];
}
}
else {
if (pflag > 1) {
fprintf(stderr," Using higher version of parameters for");
fprintf(stderr," %s ",item->keyword);
for (i=0; i < item->number_of_members; i++)
fprintf(stderr,"%s ",item->data[replace].ff_types[i]);
fprintf(stderr," version %3.2f\n",item->data[replace].ver);
}
}
}
else {
item->data[ctr].ver = version;
item->data[ctr].ref = reference;
for (i=0; i < item->number_of_members; i++) {
strncpy(item->data[ctr].ff_types[i],atom_types[i],5);
}
for (i=0; i < item->number_of_parameters; i++) {
item->data[ctr].ff_param[i] = parameters[i];
}
ctr++;
}
fgets( line, MAX_LINE, FrcF);
/*if blank line encountered, get next */
while((blank_line(line)) ||
(strncmp(line,"!",1) == 0)) {
fgets( line, MAX_LINE, FrcF);
}
}
item->entries = ctr;
/*Debugging
fprintf(stderr,"\n%s\n", item->keyword);
for(i=0;i<ctr;i++) {
for(j=0;j<item->number_of_members;j++)
fprintf(stderr,"%3s ", item->data[i].ff_equiv[j]);
fprintf(stderr," ");
for(j=0;j<item->number_of_parameters;j++)
fprintf(stderr,"%10.5f ",item->data[i].ff_param[j]);
fprintf(stderr,"\n");
}
*/
}
unsigned char string_match(char *string1,char *string2)
{
int len1,len2;
len1 = strlen(string1);
len2 = strlen(string2);
if (len1 != len2) {
return 0;
}
else {
if (strncmp(string1,string2,len1) == 0) {
return 1;
}
else {
return 0;
}
}
}

626
tools/msi2lmp/src/WriteDataFile.c
View File

@ -1,313 +1,313 @@
/*
* This function creates and writes the data file to be used with LAMMPS
*/
#include "Msi2LMP2.h"
#include "Forcefield.h"
void WriteDataFile(FILE *DatF,char *nameroot,int forcefield)
{
int i,j,k,m;
char line[MAX_LINE_LENGTH];
if (forcefield == 1) total_no_angle_angles = 0;
fprintf(DatF, "LAMMPS data file for %s\n\n", nameroot);
fprintf(DatF, " %6d atoms\n", total_no_atoms);
fprintf(DatF, " %6d bonds\n", total_no_bonds);
fprintf(DatF, " %6d angles\n",total_no_angles);
fprintf(DatF, " %6d dihedrals\n", total_no_dihedrals);
fprintf(DatF, " %6d impropers\n", total_no_oops+total_no_angle_angles);
fprintf(DatF, "\n");
fprintf(DatF, " %3d atom types\n", no_atom_types);
if (no_bond_types > 0)
fprintf(DatF, " %3d bond types\n", no_bond_types);
if (no_angle_types> 0)
fprintf(DatF, " %3d angle types\n", no_angle_types);
if (no_dihedral_types > 0) fprintf (DatF," %3d dihedral types\n",
no_dihedral_types);
if (forcefield > 1) {
if ((no_oop_types + no_angleangle_types) > 0)
fprintf (DatF, " %3d improper types\n",
no_oop_types + no_angleangle_types);
}
else {
if (no_oop_types > 0)
fprintf (DatF, " %3d improper types\n", no_oop_types);
}
fprintf(DatF, "\n");
fprintf(DatF, " %15.9f %15.9f xlo xhi\n", pbc[0], pbc[3]);
fprintf(DatF, " %15.9f %15.9f ylo yhi\n", pbc[1], pbc[4]);
fprintf(DatF, " %15.9f %15.9f zlo zhi\n", pbc[2], pbc[5]);
/* MASSES */
fprintf(DatF, "\nMasses\n\n");
for(k=0; k < no_atom_types; k++)
fprintf(DatF, " %3d %10.6f\n",k+1,atomtypes[k].mass);
fprintf(DatF, "\n");
/* COEFFICIENTS */
fprintf(DatF,"Nonbond Coeffs\n\n");
for (i=0; i < no_atom_types; i++) {
fprintf(DatF, " %3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%14.10f ", atomtypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
if (no_bond_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Bond Coeffs\n\n");
for (i=0; i < no_bond_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", bondtypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_angle_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Angle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", angletypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
if (forcefield == 1) m = 3;
if (forcefield == 2) m = 6;
fprintf(DatF,"Dihedral Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", dihedraltypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (forcefield == 1) {
if (no_oop_types > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
else if (forcefield == 2) {
if ((no_oop_types + no_angleangle_types) > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", 0.0);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
if (forcefield == 2) {
if (no_angle_types > 0) {
fprintf(DatF,"BondBond Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", angletypes[i].bondbond_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondAngle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF, "%10.4f ",angletypes[i].bondangle_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if ((no_oop_types+no_angleangle_types) > 0) {
fprintf(DatF,"AngleAngle Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].angleangle_params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", angleangletypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
fprintf(DatF,"AngleAngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].angleangledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"EndBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].endbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"MiddleBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].midbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondBond13 Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].bond13_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"AngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].angledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
/*--------------------------------------------------------------------*/
/* ATOMS */
fprintf(DatF, "Atoms\n\n");
for(k=0; k < total_no_atoms; k++) {
fprintf(DatF, " %6i %6i %3i %9.6f %15.9f %15.9f %15.9f\n",
k+1,
atoms[k].molecule,
atoms[k].type+1,
atoms[k].q,
atoms[k].x[0],
atoms[k].x[1],
atoms[k].x[2]);
}
fprintf(DatF, "\n");
/***** BONDS *****/
if (total_no_bonds > 0) {
fprintf(DatF, "Bonds\n\n");
for(k=0; k < total_no_bonds; k++)
fprintf(DatF, "%6i %3i %6i %6i\n",k+1,
bonds[k].type+1,
bonds[k].members[0]+1,
bonds[k].members[1]+1);
fprintf(DatF,"\n");
}
/***** ANGLES *****/
if (total_no_angles > 0) {
fprintf(DatF, "Angles\n\n");
for(k=0; k < total_no_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i\n",k+1,
angles[k].type+1,
angles[k].members[0]+1,
angles[k].members[1]+1,
angles[k].members[2]+1);
fprintf(DatF, "\n");
}
/***** TORSIONS *****/
if (total_no_dihedrals > 0) {
fprintf(DatF,"Dihedrals\n\n");
for(k=0; k < total_no_dihedrals; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i\n",k+1,
dihedrals[k].type+1,
dihedrals[k].members[0]+1,
dihedrals[k].members[1]+1,
dihedrals[k].members[2]+1,
dihedrals[k].members[3]+1);
fprintf(DatF, "\n");
}
/***** OUT-OF-PLANES *****/
if (total_no_oops+total_no_angle_angles > 0) {
fprintf(DatF,"Impropers\n\n");
for (k=0; k < total_no_oops; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n", k+1,
oops[k].type+1,
oops[k].members[0]+1,
oops[k].members[1]+1,
oops[k].members[2]+1,
oops[k].members[3]+1);
if (forcefield == 2) {
for (k=0; k < total_no_angle_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n",k+total_no_oops+1,
angleangles[k].type+no_oop_types+1,
angleangles[k].members[0]+1,
angleangles[k].members[1]+1,
angleangles[k].members[2]+1,
angleangles[k].members[3]+1);
}
fprintf(DatF, "\n");
}
}
/*
* This function creates and writes the data file to be used with LAMMPS
*/
#include "Msi2LMP2.h"
#include "Forcefield.h"
void WriteDataFile(FILE *DatF,char *nameroot,int forcefield)
{
int i,j,k,m;
char line[MAX_LINE_LENGTH];
if (forcefield == 1) total_no_angle_angles = 0;
fprintf(DatF, "LAMMPS data file for %s\n\n", nameroot);
fprintf(DatF, " %6d atoms\n", total_no_atoms);
fprintf(DatF, " %6d bonds\n", total_no_bonds);
fprintf(DatF, " %6d angles\n",total_no_angles);
fprintf(DatF, " %6d dihedrals\n", total_no_dihedrals);
fprintf(DatF, " %6d impropers\n", total_no_oops+total_no_angle_angles);
fprintf(DatF, "\n");
fprintf(DatF, " %3d atom types\n", no_atom_types);
if (no_bond_types > 0)
fprintf(DatF, " %3d bond types\n", no_bond_types);
if (no_angle_types> 0)
fprintf(DatF, " %3d angle types\n", no_angle_types);
if (no_dihedral_types > 0) fprintf (DatF," %3d dihedral types\n",
no_dihedral_types);
if (forcefield > 1) {
if ((no_oop_types + no_angleangle_types) > 0)
fprintf (DatF, " %3d improper types\n",
no_oop_types + no_angleangle_types);
}
else {
if (no_oop_types > 0)
fprintf (DatF, " %3d improper types\n", no_oop_types);
}
fprintf(DatF, "\n");
fprintf(DatF, " %15.9f %15.9f xlo xhi\n", pbc[0], pbc[3]);
fprintf(DatF, " %15.9f %15.9f ylo yhi\n", pbc[1], pbc[4]);
fprintf(DatF, " %15.9f %15.9f zlo zhi\n", pbc[2], pbc[5]);
/* MASSES */
fprintf(DatF, "\nMasses\n\n");
for(k=0; k < no_atom_types; k++)
fprintf(DatF, " %3d %10.6f\n",k+1,atomtypes[k].mass);
fprintf(DatF, "\n");
/* COEFFICIENTS */
fprintf(DatF,"Nonbond Coeffs\n\n");
for (i=0; i < no_atom_types; i++) {
fprintf(DatF, " %3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%14.10f ", atomtypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
if (no_bond_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Bond Coeffs\n\n");
for (i=0; i < no_bond_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", bondtypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_angle_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Angle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", angletypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
if (forcefield == 1) m = 3;
if (forcefield == 2) m = 6;
fprintf(DatF,"Dihedral Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", dihedraltypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (forcefield == 1) {
if (no_oop_types > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
else if (forcefield == 2) {
if ((no_oop_types + no_angleangle_types) > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", 0.0);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
if (forcefield == 2) {
if (no_angle_types > 0) {
fprintf(DatF,"BondBond Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", angletypes[i].bondbond_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondAngle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF, "%10.4f ",angletypes[i].bondangle_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if ((no_oop_types+no_angleangle_types) > 0) {
fprintf(DatF,"AngleAngle Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].angleangle_params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", angleangletypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
fprintf(DatF,"AngleAngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].angleangledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"EndBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].endbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"MiddleBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].midbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondBond13 Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].bond13_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"AngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].angledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
/*--------------------------------------------------------------------*/
/* ATOMS */
fprintf(DatF, "Atoms\n\n");
for(k=0; k < total_no_atoms; k++) {
fprintf(DatF, " %6i %6i %3i %9.6f %15.9f %15.9f %15.9f\n",
k+1,
atoms[k].molecule,
atoms[k].type+1,
atoms[k].q,
atoms[k].x[0],
atoms[k].x[1],
atoms[k].x[2]);
}
fprintf(DatF, "\n");
/***** BONDS *****/
if (total_no_bonds > 0) {
fprintf(DatF, "Bonds\n\n");
for(k=0; k < total_no_bonds; k++)
fprintf(DatF, "%6i %3i %6i %6i\n",k+1,
bonds[k].type+1,
bonds[k].members[0]+1,
bonds[k].members[1]+1);
fprintf(DatF,"\n");
}
/***** ANGLES *****/
if (total_no_angles > 0) {
fprintf(DatF, "Angles\n\n");
for(k=0; k < total_no_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i\n",k+1,
angles[k].type+1,
angles[k].members[0]+1,
angles[k].members[1]+1,
angles[k].members[2]+1);
fprintf(DatF, "\n");
}
/***** TORSIONS *****/
if (total_no_dihedrals > 0) {
fprintf(DatF,"Dihedrals\n\n");
for(k=0; k < total_no_dihedrals; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i\n",k+1,
dihedrals[k].type+1,
dihedrals[k].members[0]+1,
dihedrals[k].members[1]+1,
dihedrals[k].members[2]+1,
dihedrals[k].members[3]+1);
fprintf(DatF, "\n");
}
/***** OUT-OF-PLANES *****/
if (total_no_oops+total_no_angle_angles > 0) {
fprintf(DatF,"Impropers\n\n");
for (k=0; k < total_no_oops; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n", k+1,
oops[k].type+1,
oops[k].members[0]+1,
oops[k].members[1]+1,
oops[k].members[2]+1,
oops[k].members[3]+1);
if (forcefield == 2) {
for (k=0; k < total_no_angle_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n",k+total_no_oops+1,
angleangles[k].type+no_oop_types+1,
angleangles[k].members[0]+1,
angleangles[k].members[1]+1,
angleangles[k].members[2]+1,
angleangles[k].members[3]+1);
}
fprintf(DatF, "\n");
}
}

664
tools/msi2lmp/src/WriteDataFile01.c
View File

@ -1,332 +1,332 @@
/*
* This function creates and writes the data file to be used with LAMMPS
*/
#include "Msi2LMP2.h"
#include "Forcefield.h"
void WriteDataFile01(char *nameroot,int forcefield)
{
int i,j,k,m;
char line[MAX_LINE_LENGTH];
FILE *DatF;
/* Open data file */
sprintf(line,"%s.lammps01",rootname);
if (pflag > 0) fprintf(stderr," Writing LAMMPS 2001 data file: %s\n",line);
if( (DatF = fopen(line,"w")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
if (forcefield == 1) total_no_angle_angles = 0;
fprintf(DatF, "LAMMPS data file for %s\n\n", nameroot);
fprintf(DatF, " %6d atoms\n", total_no_atoms);
fprintf(DatF, " %6d bonds\n", total_no_bonds);
fprintf(DatF, " %6d angles\n",total_no_angles);
fprintf(DatF, " %6d dihedrals\n", total_no_dihedrals);
fprintf(DatF, " %6d impropers\n", total_no_oops+total_no_angle_angles);
fprintf(DatF, "\n");
fprintf(DatF, " %3d atom types\n", no_atom_types);
if (no_bond_types > 0)
fprintf(DatF, " %3d bond types\n", no_bond_types);
if (no_angle_types> 0)
fprintf(DatF, " %3d angle types\n", no_angle_types);
if (no_dihedral_types > 0) fprintf (DatF," %3d dihedral types\n",
no_dihedral_types);
if (forcefield > 1) {
if ((no_oop_types + no_angleangle_types) > 0)
fprintf (DatF, " %3d improper types\n",
no_oop_types + no_angleangle_types);
}
else {
if (no_oop_types > 0)
fprintf (DatF, " %3d improper types\n", no_oop_types);
}
fprintf(DatF, "\n");
fprintf(DatF, " %15.9f %15.9f xlo xhi\n", pbc[0], pbc[3]);
fprintf(DatF, " %15.9f %15.9f ylo yhi\n", pbc[1], pbc[4]);
fprintf(DatF, " %15.9f %15.9f zlo zhi\n", pbc[2], pbc[5]);
/* MASSES */
fprintf(DatF, "\nMasses\n\n");
for(k=0; k < no_atom_types; k++)
fprintf(DatF, " %3d %10.6f\n",k+1,atomtypes[k].mass);
fprintf(DatF, "\n");
/* COEFFICIENTS */
fprintf(DatF,"Nonbond Coeffs\n\n");
for (i=0; i < no_atom_types; i++) {
fprintf(DatF, " %3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%14.10f ", atomtypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
if (no_bond_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Bond Coeffs\n\n");
for (i=0; i < no_bond_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", bondtypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_angle_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Angle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", angletypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
if (forcefield == 1) m = 3;
if (forcefield == 2) m = 6;
fprintf(DatF,"Dihedral Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", dihedraltypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (forcefield == 1) {
if (no_oop_types > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
else if (forcefield == 2) {
if ((no_oop_types + no_angleangle_types) > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", 0.0);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
if (forcefield == 2) {
if (no_angle_types > 0) {
fprintf(DatF,"BondBond Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", angletypes[i].bondbond_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondAngle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF, "%10.4f ",angletypes[i].bondangle_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if ((no_oop_types+no_angleangle_types) > 0) {
fprintf(DatF,"AngleAngle Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].angleangle_params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", angleangletypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
fprintf(DatF,"AngleAngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].angleangledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"EndBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].endbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"MiddleBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].midbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondBond13 Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].bond13_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"AngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].angledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
/*--------------------------------------------------------------------*/
/* ATOMS */
fprintf(DatF, "Atoms\n\n");
for(k=0; k < total_no_atoms; k++) {
fprintf(DatF, " %6i %6i %3i %9.6f %15.9f %15.9f %15.9f\n",
k+1,
atoms[k].molecule,
atoms[k].type+1,
atoms[k].q,
atoms[k].x[0],
atoms[k].x[1],
atoms[k].x[2]);
}
fprintf(DatF, "\n");
/***** BONDS *****/
if (total_no_bonds > 0) {
fprintf(DatF, "Bonds\n\n");
for(k=0; k < total_no_bonds; k++)
fprintf(DatF, "%6i %3i %6i %6i\n",k+1,
bonds[k].type+1,
bonds[k].members[0]+1,
bonds[k].members[1]+1);
fprintf(DatF,"\n");
}
/***** ANGLES *****/
if (total_no_angles > 0) {
fprintf(DatF, "Angles\n\n");
for(k=0; k < total_no_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i\n",k+1,
angles[k].type+1,
angles[k].members[0]+1,
angles[k].members[1]+1,
angles[k].members[2]+1);
fprintf(DatF, "\n");
}
/***** TORSIONS *****/
if (total_no_dihedrals > 0) {
fprintf(DatF,"Dihedrals\n\n");
for(k=0; k < total_no_dihedrals; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i\n",k+1,
dihedrals[k].type+1,
dihedrals[k].members[0]+1,
dihedrals[k].members[1]+1,
dihedrals[k].members[2]+1,
dihedrals[k].members[3]+1);
fprintf(DatF, "\n");
}
/***** OUT-OF-PLANES *****/
if (total_no_oops+total_no_angle_angles > 0) {
fprintf(DatF,"Impropers\n\n");
for (k=0; k < total_no_oops; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n", k+1,
oops[k].type+1,
oops[k].members[0]+1,
oops[k].members[1]+1,
oops[k].members[2]+1,
oops[k].members[3]+1);
if (forcefield == 2) {
for (k=0; k < total_no_angle_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n",k+total_no_oops+1,
angleangles[k].type+no_oop_types+1,
angleangles[k].members[0]+1,
angleangles[k].members[1]+1,
angleangles[k].members[2]+1,
angleangles[k].members[3]+1);
}
fprintf(DatF, "\n");
}
/* Close data file */
if (fclose(DatF) !=0) {
fprintf(stderr,"Error closing %s.lammps01\n", rootname);
exit(1);
}
}
/*
* This function creates and writes the data file to be used with LAMMPS
*/
#include "Msi2LMP2.h"
#include "Forcefield.h"
void WriteDataFile01(char *nameroot,int forcefield)
{
int i,j,k,m;
char line[MAX_LINE_LENGTH];
FILE *DatF;
/* Open data file */
sprintf(line,"%s.lammps01",rootname);
if (pflag > 0) fprintf(stderr," Writing LAMMPS 2001 data file: %s\n",line);
if( (DatF = fopen(line,"w")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
if (forcefield == 1) total_no_angle_angles = 0;
fprintf(DatF, "LAMMPS data file for %s\n\n", nameroot);
fprintf(DatF, " %6d atoms\n", total_no_atoms);
fprintf(DatF, " %6d bonds\n", total_no_bonds);
fprintf(DatF, " %6d angles\n",total_no_angles);
fprintf(DatF, " %6d dihedrals\n", total_no_dihedrals);
fprintf(DatF, " %6d impropers\n", total_no_oops+total_no_angle_angles);
fprintf(DatF, "\n");
fprintf(DatF, " %3d atom types\n", no_atom_types);
if (no_bond_types > 0)
fprintf(DatF, " %3d bond types\n", no_bond_types);
if (no_angle_types> 0)
fprintf(DatF, " %3d angle types\n", no_angle_types);
if (no_dihedral_types > 0) fprintf (DatF," %3d dihedral types\n",
no_dihedral_types);
if (forcefield > 1) {
if ((no_oop_types + no_angleangle_types) > 0)
fprintf (DatF, " %3d improper types\n",
no_oop_types + no_angleangle_types);
}
else {
if (no_oop_types > 0)
fprintf (DatF, " %3d improper types\n", no_oop_types);
}
fprintf(DatF, "\n");
fprintf(DatF, " %15.9f %15.9f xlo xhi\n", pbc[0], pbc[3]);
fprintf(DatF, " %15.9f %15.9f ylo yhi\n", pbc[1], pbc[4]);
fprintf(DatF, " %15.9f %15.9f zlo zhi\n", pbc[2], pbc[5]);
/* MASSES */
fprintf(DatF, "\nMasses\n\n");
for(k=0; k < no_atom_types; k++)
fprintf(DatF, " %3d %10.6f\n",k+1,atomtypes[k].mass);
fprintf(DatF, "\n");
/* COEFFICIENTS */
fprintf(DatF,"Nonbond Coeffs\n\n");
for (i=0; i < no_atom_types; i++) {
fprintf(DatF, " %3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%14.10f ", atomtypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
if (no_bond_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Bond Coeffs\n\n");
for (i=0; i < no_bond_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", bondtypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_angle_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Angle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", angletypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
if (forcefield == 1) m = 3;
if (forcefield == 2) m = 6;
fprintf(DatF,"Dihedral Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", dihedraltypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (forcefield == 1) {
if (no_oop_types > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
else if (forcefield == 2) {
if ((no_oop_types + no_angleangle_types) > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", 0.0);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
if (forcefield == 2) {
if (no_angle_types > 0) {
fprintf(DatF,"BondBond Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", angletypes[i].bondbond_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondAngle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF, "%10.4f ",angletypes[i].bondangle_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if ((no_oop_types+no_angleangle_types) > 0) {
fprintf(DatF,"AngleAngle Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].angleangle_params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", angleangletypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
fprintf(DatF,"AngleAngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].angleangledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"EndBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].endbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"MiddleBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].midbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondBond13 Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].bond13_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"AngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].angledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
/*--------------------------------------------------------------------*/
/* ATOMS */
fprintf(DatF, "Atoms\n\n");
for(k=0; k < total_no_atoms; k++) {
fprintf(DatF, " %6i %6i %3i %9.6f %15.9f %15.9f %15.9f\n",
k+1,
atoms[k].molecule,
atoms[k].type+1,
atoms[k].q,
atoms[k].x[0],
atoms[k].x[1],
atoms[k].x[2]);
}
fprintf(DatF, "\n");
/***** BONDS *****/
if (total_no_bonds > 0) {
fprintf(DatF, "Bonds\n\n");
for(k=0; k < total_no_bonds; k++)
fprintf(DatF, "%6i %3i %6i %6i\n",k+1,
bonds[k].type+1,
bonds[k].members[0]+1,
bonds[k].members[1]+1);
fprintf(DatF,"\n");
}
/***** ANGLES *****/
if (total_no_angles > 0) {
fprintf(DatF, "Angles\n\n");
for(k=0; k < total_no_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i\n",k+1,
angles[k].type+1,
angles[k].members[0]+1,
angles[k].members[1]+1,
angles[k].members[2]+1);
fprintf(DatF, "\n");
}
/***** TORSIONS *****/
if (total_no_dihedrals > 0) {
fprintf(DatF,"Dihedrals\n\n");
for(k=0; k < total_no_dihedrals; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i\n",k+1,
dihedrals[k].type+1,
dihedrals[k].members[0]+1,
dihedrals[k].members[1]+1,
dihedrals[k].members[2]+1,
dihedrals[k].members[3]+1);
fprintf(DatF, "\n");
}
/***** OUT-OF-PLANES *****/
if (total_no_oops+total_no_angle_angles > 0) {
fprintf(DatF,"Impropers\n\n");
for (k=0; k < total_no_oops; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n", k+1,
oops[k].type+1,
oops[k].members[0]+1,
oops[k].members[1]+1,
oops[k].members[2]+1,
oops[k].members[3]+1);
if (forcefield == 2) {
for (k=0; k < total_no_angle_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n",k+total_no_oops+1,
angleangles[k].type+no_oop_types+1,
angleangles[k].members[0]+1,
angleangles[k].members[1]+1,
angleangles[k].members[2]+1,
angleangles[k].members[3]+1);
}
fprintf(DatF, "\n");
}
/* Close data file */
if (fclose(DatF) !=0) {
fprintf(stderr,"Error closing %s.lammps01\n", rootname);
exit(1);
}
}

666
tools/msi2lmp/src/WriteDataFile05.c
View File

@ -1,333 +1,333 @@
/*
* This function creates and writes the data file to be used with LAMMPS
*/
#include "Msi2LMP2.h"
#include "Forcefield.h"
void WriteDataFile05(char *nameroot,int forcefield)
{
int i,j,k,m;
char line[MAX_LINE_LENGTH];
FILE *DatF;
/* Open data file */
sprintf(line,"%s.lammps05",rootname);
if (pflag > 0) fprintf(stderr," Writing LAMMPS 2005 data file: %s\n",line);
if( (DatF = fopen(line,"w")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
if (forcefield == 1) total_no_angle_angles = 0;
fprintf(DatF, "LAMMPS 2005 data file for %s\n\n", nameroot);
fprintf(DatF, " %6d atoms\n", total_no_atoms);
fprintf(DatF, " %6d bonds\n", total_no_bonds);
fprintf(DatF, " %6d angles\n",total_no_angles);
fprintf(DatF, " %6d dihedrals\n", total_no_dihedrals);
fprintf(DatF, " %6d impropers\n", total_no_oops+total_no_angle_angles);
fprintf(DatF, "\n");
fprintf(DatF, " %3d atom types\n", no_atom_types);
if (no_bond_types > 0)
fprintf(DatF, " %3d bond types\n", no_bond_types);
if (no_angle_types> 0)
fprintf(DatF, " %3d angle types\n", no_angle_types);
if (no_dihedral_types > 0) fprintf (DatF," %3d dihedral types\n",
no_dihedral_types);
if (forcefield > 1) {
if ((no_oop_types + no_angleangle_types) > 0)
fprintf (DatF, " %3d improper types\n",
no_oop_types + no_angleangle_types);
}
else {
if (no_oop_types > 0)
fprintf (DatF, " %3d improper types\n", no_oop_types);
}
fprintf(DatF, "\n");
fprintf(DatF, " %15.9f %15.9f xlo xhi\n", pbc[0], pbc[3]);
fprintf(DatF, " %15.9f %15.9f ylo yhi\n", pbc[1], pbc[4]);
fprintf(DatF, " %15.9f %15.9f zlo zhi\n", pbc[2], pbc[5]);
/* MASSES */
fprintf(DatF, "\nMasses\n\n");
for(k=0; k < no_atom_types; k++)
fprintf(DatF, " %3d %10.6f\n",k+1,atomtypes[k].mass);
fprintf(DatF, "\n");
/* COEFFICIENTS */
fprintf(DatF,"Pair Coeffs\n\n");
for (i=0; i < no_atom_types; i++) {
fprintf(DatF, " %3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%14.10f ", atomtypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
if (no_bond_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Bond Coeffs\n\n");
for (i=0; i < no_bond_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", bondtypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_angle_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Angle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", angletypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
if (forcefield == 1) m = 3;
if (forcefield == 2) m = 6;
fprintf(DatF,"Dihedral Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", dihedraltypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (forcefield == 1) {
if (no_oop_types > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
else if (forcefield == 2) {
if ((no_oop_types + no_angleangle_types) > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", 0.0);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
if (forcefield == 2) {
if (no_angle_types > 0) {
fprintf(DatF,"BondBond Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", angletypes[i].bondbond_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondAngle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF, "%10.4f ",angletypes[i].bondangle_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if ((no_oop_types+no_angleangle_types) > 0) {
fprintf(DatF,"AngleAngle Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].angleangle_params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", angleangletypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
fprintf(DatF,"AngleAngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].angleangledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"EndBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].endbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"MiddleBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].midbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondBond13 Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].bond13_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"AngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].angledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
/*--------------------------------------------------------------------*/
/* ATOMS */
fprintf(DatF, "Atoms\n\n");
for(k=0; k < total_no_atoms; k++) {
fprintf(DatF, " %6i %6i %3i %9.6f %15.9f %15.9f %15.9f\n",
k+1,
atoms[k].molecule,
atoms[k].type+1,
atoms[k].q,
atoms[k].x[0],
atoms[k].x[1],
atoms[k].x[2]);
}
fprintf(DatF, "\n");
/***** BONDS *****/
if (total_no_bonds > 0) {
fprintf(DatF, "Bonds\n\n");
for(k=0; k < total_no_bonds; k++)
fprintf(DatF, "%6i %3i %6i %6i\n",k+1,
bonds[k].type+1,
bonds[k].members[0]+1,
bonds[k].members[1]+1);
fprintf(DatF,"\n");
}
/***** ANGLES *****/
if (total_no_angles > 0) {
fprintf(DatF, "Angles\n\n");
for(k=0; k < total_no_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i\n",k+1,
angles[k].type+1,
angles[k].members[0]+1,
angles[k].members[1]+1,
angles[k].members[2]+1);
fprintf(DatF, "\n");
}
/***** TORSIONS *****/
if (total_no_dihedrals > 0) {
fprintf(DatF,"Dihedrals\n\n");
for(k=0; k < total_no_dihedrals; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i\n",k+1,
dihedrals[k].type+1,
dihedrals[k].members[0]+1,
dihedrals[k].members[1]+1,
dihedrals[k].members[2]+1,
dihedrals[k].members[3]+1);
fprintf(DatF, "\n");
}
/***** OUT-OF-PLANES *****/
if (total_no_oops+total_no_angle_angles > 0) {
fprintf(DatF,"Impropers\n\n");
for (k=0; k < total_no_oops; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n", k+1,
oops[k].type+1,
oops[k].members[0]+1,
oops[k].members[1]+1,
oops[k].members[2]+1,
oops[k].members[3]+1);
if (forcefield == 2) {
for (k=0; k < total_no_angle_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n",k+total_no_oops+1,
angleangles[k].type+no_oop_types+1,
angleangles[k].members[0]+1,
angleangles[k].members[1]+1,
angleangles[k].members[2]+1,
angleangles[k].members[3]+1);
}
fprintf(DatF, "\n");
}
/* Close data file */
if (fclose(DatF) !=0) {
fprintf(stderr,"Error closing %s.lammps05\n", rootname);
exit(1);
}
}
/*
* This function creates and writes the data file to be used with LAMMPS
*/
#include "Msi2LMP2.h"
#include "Forcefield.h"
void WriteDataFile05(char *nameroot,int forcefield)
{
int i,j,k,m;
char line[MAX_LINE_LENGTH];
FILE *DatF;
/* Open data file */
sprintf(line,"%s.lammps05",rootname);
if (pflag > 0) fprintf(stderr," Writing LAMMPS 2005 data file: %s\n",line);
if( (DatF = fopen(line,"w")) == NULL ) {
fprintf(stderr,"Cannot open %s\n",line);
exit(2);
}
if (forcefield == 1) total_no_angle_angles = 0;
fprintf(DatF, "LAMMPS 2005 data file for %s\n\n", nameroot);
fprintf(DatF, " %6d atoms\n", total_no_atoms);
fprintf(DatF, " %6d bonds\n", total_no_bonds);
fprintf(DatF, " %6d angles\n",total_no_angles);
fprintf(DatF, " %6d dihedrals\n", total_no_dihedrals);
fprintf(DatF, " %6d impropers\n", total_no_oops+total_no_angle_angles);
fprintf(DatF, "\n");
fprintf(DatF, " %3d atom types\n", no_atom_types);
if (no_bond_types > 0)
fprintf(DatF, " %3d bond types\n", no_bond_types);
if (no_angle_types> 0)
fprintf(DatF, " %3d angle types\n", no_angle_types);
if (no_dihedral_types > 0) fprintf (DatF," %3d dihedral types\n",
no_dihedral_types);
if (forcefield > 1) {
if ((no_oop_types + no_angleangle_types) > 0)
fprintf (DatF, " %3d improper types\n",
no_oop_types + no_angleangle_types);
}
else {
if (no_oop_types > 0)
fprintf (DatF, " %3d improper types\n", no_oop_types);
}
fprintf(DatF, "\n");
fprintf(DatF, " %15.9f %15.9f xlo xhi\n", pbc[0], pbc[3]);
fprintf(DatF, " %15.9f %15.9f ylo yhi\n", pbc[1], pbc[4]);
fprintf(DatF, " %15.9f %15.9f zlo zhi\n", pbc[2], pbc[5]);
/* MASSES */
fprintf(DatF, "\nMasses\n\n");
for(k=0; k < no_atom_types; k++)
fprintf(DatF, " %3d %10.6f\n",k+1,atomtypes[k].mass);
fprintf(DatF, "\n");
/* COEFFICIENTS */
fprintf(DatF,"Pair Coeffs\n\n");
for (i=0; i < no_atom_types; i++) {
fprintf(DatF, " %3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%14.10f ", atomtypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
if (no_bond_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Bond Coeffs\n\n");
for (i=0; i < no_bond_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", bondtypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_angle_types > 0) {
if (forcefield == 1) m = 2;
if (forcefield == 2) m = 4;
fprintf(DatF,"Angle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", angletypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
if (forcefield == 1) m = 3;
if (forcefield == 2) m = 6;
fprintf(DatF,"Dihedral Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < m; j++)
fprintf(DatF, "%10.4f ", dihedraltypes[i].params[j]);
fprintf(DatF,"\n");
}
fprintf(DatF, "\n");
}
if (forcefield == 1) {
if (no_oop_types > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
else if (forcefield == 2) {
if ((no_oop_types + no_angleangle_types) > 0) {
fprintf(DatF,"Improper Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 2; j++)
fprintf(DatF, "%10.4f ", 0.0);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
if (forcefield == 2) {
if (no_angle_types > 0) {
fprintf(DatF,"BondBond Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ", angletypes[i].bondbond_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondAngle Coeffs\n\n");
for (i=0; i < no_angle_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF, "%10.4f ",angletypes[i].bondangle_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if ((no_oop_types+no_angleangle_types) > 0) {
fprintf(DatF,"AngleAngle Coeffs\n\n");
for (i=0; i < no_oop_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", ooptypes[i].angleangle_params[j]);
fprintf(DatF, "\n");
}
for (i=0; i < no_angleangle_types; i++) {
fprintf(DatF, "%3i ", i+no_oop_types+1);
for ( j = 0; j < 6; j++)
fprintf(DatF, "%10.4f ", angleangletypes[i].params[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
if (no_dihedral_types > 0) {
fprintf(DatF,"AngleAngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].angleangledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"EndBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].endbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"MiddleBondTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 4; j++)
fprintf(DatF,"%10.4f ",
dihedraltypes[i].midbonddihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"BondBond13 Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 3; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].bond13_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
fprintf(DatF,"AngleTorsion Coeffs\n\n");
for (i=0; i < no_dihedral_types; i++) {
fprintf(DatF, "%3i ", i+1);
for ( j = 0; j < 8; j++)
fprintf(DatF, "%10.4f ",
dihedraltypes[i].angledihedral_cross_term[j]);
fprintf(DatF, "\n");
}
fprintf(DatF, "\n");
}
}
/*--------------------------------------------------------------------*/
/* ATOMS */
fprintf(DatF, "Atoms\n\n");
for(k=0; k < total_no_atoms; k++) {
fprintf(DatF, " %6i %6i %3i %9.6f %15.9f %15.9f %15.9f\n",
k+1,
atoms[k].molecule,
atoms[k].type+1,
atoms[k].q,
atoms[k].x[0],
atoms[k].x[1],
atoms[k].x[2]);
}
fprintf(DatF, "\n");
/***** BONDS *****/
if (total_no_bonds > 0) {
fprintf(DatF, "Bonds\n\n");
for(k=0; k < total_no_bonds; k++)
fprintf(DatF, "%6i %3i %6i %6i\n",k+1,
bonds[k].type+1,
bonds[k].members[0]+1,
bonds[k].members[1]+1);
fprintf(DatF,"\n");
}
/***** ANGLES *****/
if (total_no_angles > 0) {
fprintf(DatF, "Angles\n\n");
for(k=0; k < total_no_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i\n",k+1,
angles[k].type+1,
angles[k].members[0]+1,
angles[k].members[1]+1,
angles[k].members[2]+1);
fprintf(DatF, "\n");
}
/***** TORSIONS *****/
if (total_no_dihedrals > 0) {
fprintf(DatF,"Dihedrals\n\n");
for(k=0; k < total_no_dihedrals; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i\n",k+1,
dihedrals[k].type+1,
dihedrals[k].members[0]+1,
dihedrals[k].members[1]+1,
dihedrals[k].members[2]+1,
dihedrals[k].members[3]+1);
fprintf(DatF, "\n");
}
/***** OUT-OF-PLANES *****/
if (total_no_oops+total_no_angle_angles > 0) {
fprintf(DatF,"Impropers\n\n");
for (k=0; k < total_no_oops; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n", k+1,
oops[k].type+1,
oops[k].members[0]+1,
oops[k].members[1]+1,
oops[k].members[2]+1,
oops[k].members[3]+1);
if (forcefield == 2) {
for (k=0; k < total_no_angle_angles; k++)
fprintf(DatF, "%6i %3i %6i %6i %6i %6i \n",k+total_no_oops+1,
angleangles[k].type+no_oop_types+1,
angleangles[k].members[0]+1,
angleangles[k].members[1]+1,
angleangles[k].members[2]+1,
angleangles[k].members[3]+1);
}
fprintf(DatF, "\n");
}
/* Close data file */
if (fclose(DatF) !=0) {
fprintf(stderr,"Error closing %s.lammps05\n", rootname);
exit(1);
}
}

550
tools/msi2lmp/src/msi2lmp.c
View File

@ -1,275 +1,275 @@
/*
*
* msi2lmp.exe V3.6
*
* v3.6 KLA - Changes to output to either lammps 2001 (F90 version) or to
* lammps 2005 (C++ version)
*
* v3.4 JEC - a number of minor changes due to way newline and EOF are generated
* on Materials Studio generated .car and .mdf files as well as odd
* behavior out of newer Linux IO libraries. ReadMdfFile was restructured
* in the process.
*
* v3.1 JEC - changed IO interface to standard in/out, forcefield file
* location can be indicated by environmental variable; added
* printing options, consistency checks and forcefield
* parameter versions sensitivity (highest one used)
*
* v3.0 JEC - program substantially rewritten to reduce execution time
* and be 98 % dynamic in memory use (still fixed limits on
* number of parameter types for different internal coordinate
* sets)
*
* v2.0 MDP - got internal coordinate information from mdf file and
* forcefield parameters from frc file thus eliminating
* need for Discover
*
* V1.0 SL - original version. Used .car file and internal coordinate
* information from Discover to produce LAMMPS data file.
*
* This program uses the .car and .mdf files from MSI/Biosyms's INSIGHT
* program to produce a LAMMPS data file.
*
* The program is started by supplying information at the command prompt
* according to the usage described below.
*
* USAGE: msi2lmp3 ROOTNAME {-print #} {-class #} {-frc FRC_FILE} -2001
*
* -- msi2lmp3 is the name of the executable
* -- ROOTNAME is the base name of the .car and .mdf files
*
* -- -print
* # is the print level: 0 - silent except for errors
* 1 - minimal (default)
* 2 - more verbose
* -- -class
* # is the class of forcefield to use (I = Class I e.g., CVFF)
* (II = Class II e.g., CFFx )
* default is -class I
*
* -- -frc - specifies name of the forcefield file (e.g., cff91)
*
* If the name includes a hard wired directory (i.e., if the name
* starts with . or /), then the name is used alone. Otherwise,
* the program looks for the forcefield file in $BIOSYM_LIBRARY.
* If $BIOSYM_LIBRARY is not set, then the current directory is
* used.
*
* If the file name does not include a dot after the first
* character, then .frc is appended to the name.
*
* For example, -frc cvff (assumes cvff.frc is in $BIOSYM_LIBRARY
* or .)
*
* -frc cff/cff91 (assumes cff91.frc is in
* $BIOSYM_LIBRARY/cff or ./cff)
*
* -frc /usr/local/biosym/forcefields/cff95 (absolute
* location)
*
* By default, the program uses $BIOSYM_LIBRARY/cvff.frc
*
* -- -2001 will output a data file for the FORTRAN 90 version of LAMMPS (2001)
* By default, the program will output for the C++ version of LAMMPS.
*
* -- output is written to a file called ROOTNAME.lammps{01/05}
*
*
****************************************************************
*
* Msi2lmp3
*
* This is the third version of a program that generates a LAMMPS
* data file based on the information in a MSI car file (atom
* coordinates) and mdf file (molecular topology). A key part of
* the program looks up forcefield parameters from an MSI frc file.
*
* The first version was written by Steve Lustig at Dupont, but
* required using Discover to derive internal coordinates and
* forcefield parameters
*
* The second version was written by Michael Peachey while an
* in intern in the Cray Chemistry Applications Group managed
* by John Carpenter. This version derived internal coordinates
* from the mdf file and looked up parameters in the frc file
* thus eliminating the need for Discover.
*
* The third version was written by John Carpenter to optimize
* the performance of the program for large molecular systems
* (the original code for deriving atom numbers was quadratic in time)
* and to make the program fully dynamic. The second version used
* fixed dimension arrays for the internal coordinates.
*
* John Carpenter can be contacted by sending email to
* jec374@earthlink.net
*
* November 2000
*/
#define MAIN
#include "Msi2LMP2.h"
int main (int argc, char *argv[])
{
int n,i,found_dot; /* Counter */
int outv;
char *string;
char *frc_dir_name;
char *frc_file_name;
FILE *DatF;
/* Functions called from within main */
/* All code is located in .c file with function name */
extern void FrcMenu();
extern void ReadCarFile();
extern void ReadMdfFile();
extern void ReadFrcFile();
extern void MakeLists();
extern void GetParameters(int);
extern void CheckLists();
extern void WriteDataFile(FILE *,char *,int);
outv = 2005;
pflag = 1;
forcefield = 1; /* Variable that identifies forcefield to use */
frc_file_name = (char *) calloc(160,sizeof(char));
frc_dir_name = (char *) calloc(160,sizeof(char));
frc_dir_name = getenv("BIOSYM_LIBRARY");
if (frc_dir_name == NULL) {
frc_file_name = strcpy(frc_file_name,"./cvff.frc");
}
else {
for (i=0; i < strlen(frc_dir_name); i++)
frc_file_name[i] = frc_dir_name[i];
frc_file_name = strcat(frc_file_name,"/cvff.frc");
}
if (argc < 2) { /* If no rootname was supplied, prompt for it */
fprintf(stderr,"The rootname of the .car and .mdf files must be entered\n");
}
else /* rootname was supplied as first argument, copy to rootname */
sprintf(rootname,"%s",argv[1]);
n = 2;
while (n < argc) {
if (strcmp(argv[n],"-class") == 0) {
if (strcmp(argv[n+1],"I") == 0) {
forcefield = 1;
n++;
}
else if (strcmp(argv[n+1],"II") == 0) {
forcefield = 2;
n++;
}
else {
fprintf(stderr,"Unrecognized Forcefield class: %s\n",
argv[n+1]);
n++;
}
}
else if (strcmp(argv[n],"-2001") == 0) {
outv = 2001;
n++;
}
else if (strcmp(argv[n],"-frc") == 0) {
if ((frc_dir_name != NULL) && (argv[n+1][0] != '.')) {
for (i=0; i < strlen(frc_dir_name); i++) {
frc_file_name[i] = frc_dir_name[i];
}
frc_file_name[strlen(frc_dir_name)] = '\0';
frc_file_name = strcat(frc_file_name,"/");
frc_file_name = strcat(frc_file_name,argv[n+1]);
}
else {
frc_file_name = strcpy(frc_file_name,argv[n+1]);
}
found_dot = 0;
for (i=1; i < strlen(frc_file_name); i++) {
if (frc_file_name[i] == '.') found_dot = 1;
}
if (found_dot == 0)
frc_file_name = strcat(frc_file_name,".frc");
n++;
}
else if (strstr(argv[n],"-p") != NULL) {
pflag = atoi(argv[n+1]);
n++;
}
else {
fprintf(stderr,"Unrecognized option: %s\n",argv[n]);
}
n++;
}
for (i=0; i < strlen(frc_file_name); i++)
FrcFileName[i] = frc_file_name[i];
free(frc_file_name);
if (pflag > 0) {
fprintf(stderr,"\nRunning Msi2lmp.....\n\n");
fprintf(stderr," Forcefield file name: %s\n",FrcFileName);
fprintf(stderr," Forcefield class: %d\n\n",forcefield);
}
if (((forcefield == 1) && (strstr(FrcFileName,"cff") != NULL) ||
(forcefield == 2) && (strstr(FrcFileName,"cvff") != NULL))) {
fprintf(stderr," WARNING - forcefield name and class appear to\n");
fprintf(stderr," be inconsistent - Errors may result\n\n");
}
/* Read in .car file */
ReadCarFile();
/*Read in .mdf file */
ReadMdfFile();
/* Define bonds, angles, etc...*/
if (pflag > 0) fprintf(stderr,"\n Building internal coordinate lists \n");
MakeLists();
/* Read .frc file into memory */
if (pflag > 0) fprintf(stderr,"\n Reading forcefield file \n");
ReadFrcFile();
/* Get forcefield parameters */
if (pflag > 0) fprintf(stderr,"\n Get parameters for this molecular system\n");
GetParameters(forcefield);
/* Do internal check of internal coordinate lists */
if (pflag > 0) fprintf(stderr,"\n Check parameters for internal consistency\n");
CheckLists();
if (outv == 2001) { fprintf(stderr,"\n Writing LAMMPS 2001 data file\n");
WriteDataFile01(rootname,forcefield);
}
else if (outv == 2005) {fprintf(stderr,"\n Writing LAMMPS 2005 data file\n");
WriteDataFile05(rootname,forcefield);
}
if (pflag > 0) fprintf(stderr,"\nNormal program termination\n");
}
#include <ctype.h>
int blank_line(char *line)
{
int i,n;
for (i=0,n=0; i < strlen(line); i++) {
if (isalnum((unsigned char)line[i])) n++;
}
if (n > 0) {
return(0);
}
else {
return(1);
}
}
/*
*
* msi2lmp.exe V3.6
*
* v3.6 KLA - Changes to output to either lammps 2001 (F90 version) or to
* lammps 2005 (C++ version)
*
* v3.4 JEC - a number of minor changes due to way newline and EOF are generated
* on Materials Studio generated .car and .mdf files as well as odd
* behavior out of newer Linux IO libraries. ReadMdfFile was restructured
* in the process.
*
* v3.1 JEC - changed IO interface to standard in/out, forcefield file
* location can be indicated by environmental variable; added
* printing options, consistency checks and forcefield
* parameter versions sensitivity (highest one used)
*
* v3.0 JEC - program substantially rewritten to reduce execution time
* and be 98 % dynamic in memory use (still fixed limits on
* number of parameter types for different internal coordinate
* sets)
*
* v2.0 MDP - got internal coordinate information from mdf file and
* forcefield parameters from frc file thus eliminating
* need for Discover
*
* V1.0 SL - original version. Used .car file and internal coordinate
* information from Discover to produce LAMMPS data file.
*
* This program uses the .car and .mdf files from MSI/Biosyms's INSIGHT
* program to produce a LAMMPS data file.
*
* The program is started by supplying information at the command prompt
* according to the usage described below.
*
* USAGE: msi2lmp3 ROOTNAME {-print #} {-class #} {-frc FRC_FILE} -2001
*
* -- msi2lmp3 is the name of the executable
* -- ROOTNAME is the base name of the .car and .mdf files
*
* -- -print
* # is the print level: 0 - silent except for errors
* 1 - minimal (default)
* 2 - more verbose
* -- -class
* # is the class of forcefield to use (I = Class I e.g., CVFF)
* (II = Class II e.g., CFFx )
* default is -class I
*
* -- -frc - specifies name of the forcefield file (e.g., cff91)
*
* If the name includes a hard wired directory (i.e., if the name
* starts with . or /), then the name is used alone. Otherwise,
* the program looks for the forcefield file in $BIOSYM_LIBRARY.
* If $BIOSYM_LIBRARY is not set, then the current directory is
* used.
*
* If the file name does not include a dot after the first
* character, then .frc is appended to the name.
*
* For example, -frc cvff (assumes cvff.frc is in $BIOSYM_LIBRARY
* or .)
*
* -frc cff/cff91 (assumes cff91.frc is in
* $BIOSYM_LIBRARY/cff or ./cff)
*
* -frc /usr/local/biosym/forcefields/cff95 (absolute
* location)
*
* By default, the program uses $BIOSYM_LIBRARY/cvff.frc
*
* -- -2001 will output a data file for the FORTRAN 90 version of LAMMPS (2001)
* By default, the program will output for the C++ version of LAMMPS.
*
* -- output is written to a file called ROOTNAME.lammps{01/05}
*
*
****************************************************************
*
* Msi2lmp3
*
* This is the third version of a program that generates a LAMMPS
* data file based on the information in a MSI car file (atom
* coordinates) and mdf file (molecular topology). A key part of
* the program looks up forcefield parameters from an MSI frc file.
*
* The first version was written by Steve Lustig at Dupont, but
* required using Discover to derive internal coordinates and
* forcefield parameters
*
* The second version was written by Michael Peachey while an
* in intern in the Cray Chemistry Applications Group managed
* by John Carpenter. This version derived internal coordinates
* from the mdf file and looked up parameters in the frc file
* thus eliminating the need for Discover.
*
* The third version was written by John Carpenter to optimize
* the performance of the program for large molecular systems
* (the original code for deriving atom numbers was quadratic in time)
* and to make the program fully dynamic. The second version used
* fixed dimension arrays for the internal coordinates.
*
* John Carpenter can be contacted by sending email to
* jec374@earthlink.net
*
* November 2000
*/
#define MAIN
#include "Msi2LMP2.h"
int main (int argc, char *argv[])
{
int n,i,found_dot; /* Counter */
int outv;
char *string;
char *frc_dir_name;
char *frc_file_name;
FILE *DatF;
/* Functions called from within main */
/* All code is located in .c file with function name */
extern void FrcMenu();
extern void ReadCarFile();
extern void ReadMdfFile();
extern void ReadFrcFile();
extern void MakeLists();
extern void GetParameters(int);
extern void CheckLists();
extern void WriteDataFile(FILE *,char *,int);
outv = 2005;
pflag = 1;
forcefield = 1; /* Variable that identifies forcefield to use */
frc_file_name = (char *) calloc(160,sizeof(char));
frc_dir_name = (char *) calloc(160,sizeof(char));
frc_dir_name = getenv("BIOSYM_LIBRARY");
if (frc_dir_name == NULL) {
frc_file_name = strcpy(frc_file_name,"./cvff.frc");
}
else {
for (i=0; i < strlen(frc_dir_name); i++)
frc_file_name[i] = frc_dir_name[i];
frc_file_name = strcat(frc_file_name,"/cvff.frc");
}
if (argc < 2) { /* If no rootname was supplied, prompt for it */
fprintf(stderr,"The rootname of the .car and .mdf files must be entered\n");
}
else /* rootname was supplied as first argument, copy to rootname */
sprintf(rootname,"%s",argv[1]);
n = 2;
while (n < argc) {
if (strcmp(argv[n],"-class") == 0) {
if (strcmp(argv[n+1],"I") == 0) {
forcefield = 1;
n++;
}
else if (strcmp(argv[n+1],"II") == 0) {
forcefield = 2;
n++;
}
else {
fprintf(stderr,"Unrecognized Forcefield class: %s\n",
argv[n+1]);
n++;
}
}
else if (strcmp(argv[n],"-2001") == 0) {
outv = 2001;
n++;
}
else if (strcmp(argv[n],"-frc") == 0) {
if ((frc_dir_name != NULL) && (argv[n+1][0] != '.')) {
for (i=0; i < strlen(frc_dir_name); i++) {
frc_file_name[i] = frc_dir_name[i];
}
frc_file_name[strlen(frc_dir_name)] = '\0';
frc_file_name = strcat(frc_file_name,"/");
frc_file_name = strcat(frc_file_name,argv[n+1]);
}
else {
frc_file_name = strcpy(frc_file_name,argv[n+1]);
}
found_dot = 0;
for (i=1; i < strlen(frc_file_name); i++) {
if (frc_file_name[i] == '.') found_dot = 1;
}
if (found_dot == 0)
frc_file_name = strcat(frc_file_name,".frc");
n++;
}
else if (strstr(argv[n],"-p") != NULL) {
pflag = atoi(argv[n+1]);
n++;
}
else {
fprintf(stderr,"Unrecognized option: %s\n",argv[n]);
}
n++;
}
for (i=0; i < strlen(frc_file_name); i++)
FrcFileName[i] = frc_file_name[i];
free(frc_file_name);
if (pflag > 0) {
fprintf(stderr,"\nRunning Msi2lmp.....\n\n");
fprintf(stderr," Forcefield file name: %s\n",FrcFileName);
fprintf(stderr," Forcefield class: %d\n\n",forcefield);
}
if (((forcefield == 1) && (strstr(FrcFileName,"cff") != NULL) ||
(forcefield == 2) && (strstr(FrcFileName,"cvff") != NULL))) {
fprintf(stderr," WARNING - forcefield name and class appear to\n");
fprintf(stderr," be inconsistent - Errors may result\n\n");
}
/* Read in .car file */
ReadCarFile();
/*Read in .mdf file */
ReadMdfFile();
/* Define bonds, angles, etc...*/
if (pflag > 0) fprintf(stderr,"\n Building internal coordinate lists \n");
MakeLists();
/* Read .frc file into memory */
if (pflag > 0) fprintf(stderr,"\n Reading forcefield file \n");
ReadFrcFile();
/* Get forcefield parameters */
if (pflag > 0) fprintf(stderr,"\n Get parameters for this molecular system\n");
GetParameters(forcefield);
/* Do internal check of internal coordinate lists */
if (pflag > 0) fprintf(stderr,"\n Check parameters for internal consistency\n");
CheckLists();
if (outv == 2001) { fprintf(stderr,"\n Writing LAMMPS 2001 data file\n");
WriteDataFile01(rootname,forcefield);
}
else if (outv == 2005) {fprintf(stderr,"\n Writing LAMMPS 2005 data file\n");
WriteDataFile05(rootname,forcefield);
}
if (pflag > 0) fprintf(stderr,"\nNormal program termination\n");
}
#include <ctype.h>
int blank_line(char *line)
{
int i,n;
for (i=0,n=0; i < strlen(line); i++) {
if (isalnum((unsigned char)line[i])) n++;
}
if (n > 0) {
return(0);
}
else {
return(1);
}
}