Sunday, May 31, 2015

Code Sample 1

Directory and Intermediate Image Maker


Presented here are my first fruits. This set of code was written using Notepad++ and compiled/run in MinGW. This code has the following abilities:
  • Generate a user specified number of intermediate images (directories)
  • Generate a relevant set of atomic coordinates for each image in the form of a Gaussian submission file. This is done using a linear interpolation at the moment so it's nothing all that great in my opinion.
The Code: (to look at detail I recommend copying this into Notepad++ and setting the language option to C.)


#include <sys/types.h>
#include <sys/stat.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>


//This program is meant to demonstrate the creation of a user specified number of directories
//with a file being written in each directory. Ultimately, this program will form the base of
//the directory and file creation in the NEB program, allowing each image to be contained in
//its own directory.


int main(void)
{
    //Here the variables for the program are declared. They all function as either counters
    //or switches. Several file pointers are declared as well, these are used later one for
    //keeping track of which files are in use.

   
    int i; //'i' is for image, this is the image counter
    int k;
    FILE* fp;
    FILE *fp1 = NULL;
    FILE *fp2 = NULL;
    int j;
    int number_of_atoms;
    float temp; //this is used to transfer the integer value of 'i' to a floating point number
    int c;
    c = 1;
    int q = 0;
   
    //Prompt the user for how many images they would like the NEB algorithm to use. This then
    //creates the desired number of directories.

    printf("Please enter the number of images you would like use: ");
    scanf("%d", &i);

    //set temp equal to i. This will become important later on when determining the slope
    //for the linear interpolation between points.
    temp = i;

   
    //Open the REACTANT file.
    FILE* myfile = fopen("REACTANT.txt", "r");
    //Another counter, this one used for helping to determine the number of lines in the REACTANT.
    int ch;
    int number_of_lines = 0;
   
    do
    {
        ch = fgetc(myfile);
        if(ch == '\n')
            number_of_lines++;
    } while (ch != EOF);
   
     //last line doesn't end with a new line!
     //but there has to be a line at least before the last line

    if(ch != '\n' && number_of_lines != 0)
        number_of_lines++;
   
    fclose(myfile);
   
    printf("number of lines in REACTANT.txt = %d \n", number_of_lines);
   
    //The number of atoms in the REACTANT file are calculated here.
    number_of_atoms = number_of_lines - 6;
   
    //The position arrays for the atoms are declared here since declaring them earlier wouldn't
    //have taken into account the unknown size of the file.

    float xyz1[number_of_atoms][3];
    float xyz2[number_of_atoms][3];
    float slope[number_of_atoms][3];
    //The Interpolation Array has to be big enough to store 3 coordinate for each atom in each image.
    //Hence the calculation and its inevitable large size.

    float InterpArray[number_of_atoms*i][3];
   
    //These are the input strings used to hold the information in both REACTANT and  PRODUCT
    //files as a series of strings.
    char inputString1[(4*number_of_atoms)+10][100];
    char inputString2[(4*number_of_atoms)+10][100];
   
    //This for loop reads the REACTANT file into the first input string.
    fp1 = fopen("REACTANT.txt", "r");
    for(i = 1; i < ((4*number_of_atoms)+10); i++)
    {
        fscanf(fp1, "%s", &inputString1[i]);
    }
   
    //This for loop does so for the PRODUCT file and the second input string.
    fp2 = fopen("PRODUCT.txt", "r");
    for(i = 1; i < ((4*number_of_atoms)+10); i++)
    {
        fscanf(fp2, "%s", &inputString2[i]);
    }
   
    //Reset counter values
    k = 0;
    i=0;
    //This fills in the xzy arrays with the atomic coordinates from the input strings.
    while(k < number_of_atoms)

    {
        //IMPORTANT!!! The 9 here has its value because that is where the coordinates start on
        //the input string. If the number of arguments change in the input file, this number
        //and the numbers below it must change as well. This will be changed in the final
        //code since it is assumed that the NEB algorithm will only require 'sp' calculations.

        xyz1[k][0] = atof(inputString1[i+9]);
        xyz1[k][1] = atof(inputString1[i+10]);
        xyz1[k][2] = atof(inputString1[i+11]);
       
        xyz2[k][0] = atof(inputString2[i+9]);
        xyz2[k][1] = atof(inputString2[i+10]);
        xyz2[k][2] = atof(inputString2[i+11]);
       
        i+=4;
        k++;
    }
    k = 0;
    //Calculate the slope for each atom, this is the linear interpolation bit.
    while(k < number_of_atoms)
    {
        slope[k][0] = (xyz2[k][0] - xyz1[k][0])/(temp+1);
        slope[k][1] = (xyz2[k][1] - xyz1[k][1])/(temp+1);
        slope[k][2] = (xyz2[k][2] - xyz1[k][2])/(temp+1);
        k++;
    }   
   
    i = temp;
    int g;
    //Calculate the intermediate positions for all of the images and their corresponding atoms.
    while(q < i*number_of_atoms)
    {
        g = 0;
        while(g < number_of_atoms)
        {
            InterpArray[q+g][0] = ((c*slope[g][0]) + xyz1[g][0]);
            InterpArray[q+g][1] = ((c*slope[g][1]) + xyz1[g][1]);
            InterpArray[q+g][2] = ((c*slope[g][2]) + xyz1[g][2]);
            g++;
        }   
        c++;
        q+=number_of_atoms;
    }   
   
    //Create the character array which will be used to store the names of the needed directories.
    //In the future it might be a good idea to take the floor/log of the number of images to calculate
    //the necessary length of each row. This would mean that the program could use an unlimited
    //number of images.

    char ImageArray[i][50];
    char WrittingArray[7][100];
    c = 0;
   
    //Here, the directories are named and stored in the array 'ImageArray.'
    while(c < i)
    {
        sprintf(ImageArray[c], "Image_%d", c + 1);
        c++;
    }

    //Reset the value of c. This is so it can be used in the next while loop.
    c = 0;
   
    //The directories are created here. 'i' directories are created as specified by the user.
    //They are produced using the 'mkdir' command which then passes its result value to 'j.'
    //This value is used to flag if the command executed successfully or not.

    while(c < i)
    {
        //Invocation of mkdir, with appropriate directory name being passed to it.
        j = mkdir(ImageArray[c]);
   
        //Check if directory creation was successful.
        if(j != 0)
        {
            printf("Directory Creation Failed For Image: %d", c);
        }
       
        //Add to c, moving to the next directory to be created.
        c++;
    }

    //Reset value of c to zero.
    c = 0;
   
    //This section of code switches between the newly created directories and adds a file to each
    //one.

    j = 0;
    while(c < i)
    {
        //Here, the program switches its working directory. This is done in the while loop, switching
        //to a new directory until the end of the images is reached.

        k = chdir(ImageArray[c]);
        //Check that the directory switch occurred.
        if(k != 0)
        {
            printf("Directory Change Failed!");
        }
        //Create the test file in the new directory.
        else
        {
            fp = fopen("TestFile.txt", "w");
            if(!fp)
            {
                printf("File Creation Failed");
            }
           
            //Here we write into the Writing Array and subsequently write the information to the
            //new file within the image directory.

            sprintf(WrittingArray[0], "%s %s %s %s \n", inputString1[1], inputString1[2],
            inputString1[3], inputString1[4]);
            sprintf(WrittingArray[1], "\n");
            sprintf(WrittingArray[2], "%s \n", inputString1[5]);
            sprintf(WrittingArray[3], "\n");
            sprintf(WrittingArray[4], "%s %s \n", inputString1[6], inputString1[7]);
            sprintf(WrittingArray[5], "\n");
            int h;
            for(h = 0; h < 5; h++)
            {
                fprintf(fp, WrittingArray[h]);
            }   
           
            k = 0;
            int m;
            m = 0;
            //Write all the atomic coordinates to the file.
            while(k < number_of_atoms)
            {
                sprintf(WrittingArray[0], "%s    %f    %f    %f \n", inputString1[8+m],
                InterpArray[j+k][0], InterpArray[j+k][1], InterpArray[j+k][2]);
                fprintf(fp, WrittingArray[0]);
                k++;
                m+=4;
            }
        }
        //Change directories to the directory above current directory.
        k = chdir("../");
        if(k != 0)
        {
            printf("Directory Change Failed");
        }
        //Increment c to switch to the next directory.
        c++;
        j+=number_of_atoms;
    }
    //Reset all counters
    j = 0;
    g = 0;
    c = 0;
    //Print off the Images and the atomic coordinates (minus elements) to the screen.
    while(j < i*number_of_atoms)
    {
        printf("Image %d: \n", (c+1));
        g = 0;
        while(g < number_of_atoms)
        {
            printf("%f %f %f \n", InterpArray[j+g][0], InterpArray[j+g][1], InterpArray[j+g][2]);
            g++;
        }   
        j += number_of_atoms;
        c++;
    }   

    return 0;   
}

Planned additions to this code:
  • Produce an XYZ file which can be used to generate an "inspection" GIF for Jmol or some other program
  • The NEB algorithm (perhaps as a separate function?)
  • Other features (they will be added as they are needed or recommended)

Saturday, May 30, 2015

Code Progress Update 1

Accomplishments:
This week I made significant progress towards completing the initial parts of the program. Most of this work has been focused on generating the necessary files in the correct format for Gaussian 09 to use later on. To date the program has the following features:
  • Generate 'i' images between reactant and product
  • Use a formatted Gaussian 09 input file containing any number of atoms (limitation being it must be formatted in a specific way)
  • Generate relevant coordinates between reactant and product
  • Storage of all coordinates in one super-array which will be used later on for spring force calculations
 Examples:
Here is a Gaussian 09 input file, you will note that the same coordinates are repeated throughout. This is because I wanted to make a simple test file, of course in an actual calculation the atoms will not be located on the same spatial coordinate.

The Reactant Coordinates:
 #T RHF/6-31G(d) Opt Freq

Title

0 1
O    0.00000    0.00000    0.00000
H    0.00000    0.00000    0.00000
H    0.00000    0.00000    0.00000
H    1.00000    1.00000    1.00000
O    2.36581    2.36581    2.36581

The Product Coordinates:
 #T RHF/6-31G(d) Opt Freq

Title

0 1
O    5.00000    5.00000    5.00000
H    5.00000    5.00000    5.00000
H    5.00000    5.00000    5.00000
H    10.0000    10.0000    10.0000
O    12.5679    12.5679    12.5679

With these two files entered, the reactant being 'REACTANT.txt' and the product being 'PRODUCT.txt' in the same directory, the user can now enter how many intermediate images they would like the program to generate.

For this example I have selected 10 intermediate images to be generated. Now, obviously I won't show each and every text file since that would be tedious and a waste of time. Rather, the program at the moment spits out the image and the relevant coordinates (minus the element) to the terminal and writes the whole output to the file under the image's relevant directory. So, the output of the above files after selecting 10 images is like this for the terminal:

Image 1:
0.454545 0.454545 0.454545
0.454545 0.454545 0.454545
0.454545 0.454545 0.454545
1.818182 1.818182 1.818182
3.293272 3.293272 3.293272
Image 2:
0.909091 0.909091 0.909091
0.909091 0.909091 0.909091
0.909091 0.909091 0.909091
2.636364 2.636364 2.636364
4.220736 4.220736 4.220736
Image 3:
1.363636 1.363636 1.363636
1.363636 1.363636 1.363636
1.363636 1.363636 1.363636
3.454545 3.454545 3.454545
5.148198 5.148198 5.148198
... (more the same) ...
Image 10:
4.545455 4.545455 4.545455
4.545455 4.545455 4.545455
4.545455 4.545455 4.545455
9.181818 9.181818 9.181818
11.640437 11.640437 11.640437

And for image 3 looks like this in the text file:

#T RHF/6-31G(d) Opt Freq

Title

0 1
O    1.363636    1.363636    1.363636
H    1.363636    1.363636    1.363636
H    1.363636    1.363636    1.363636
H    3.454545    3.454545    3.454545
O    5.148198    5.148198    5.148198

In short order then, the program works! Well, that is to say that it can successfully generate 'i' directories with a Gaussian 09 submission file written to each directory.

Next Steps:
The goals for this coming week are to complete the following:
  • Throughly comment the code and make it more legiable
  • Begin working on the NEB function
  • Write a page explaining how NEB works in theory, both in general and how it will be specifically implemented in this program
In any case, that pretty much covers all that was accomplished since Wednesday of this last week. It is my plan to write another update on either Tuesday or Wednesday next week. Till then!

Tuesday, May 26, 2015

Introductions

Hello, I am John Eric Tiessen (Jet) and welcome to the first blog post for this development blog. To get started, I want to go over the purpose of this blog, how I will be using it and why you should be interested in it.

Purpose:
The reason for creating this blog is to help document the creation of several molecular dynamic simulation tools which I will developing over the course of this summer. These tools will be an important part of my senior research project at Valparaiso University and will form the basis of any research that I intend to carry out next year.

Blog Use:
My intention with this blog is to create an open and easily accessible platform for people to both view, criticize, inquire and, perhaps, help improve the code which I will be working on. Overall, the point is to make this work as public as possible so as nothing is left undocumented and unexplained in the development process.

The Interest:
If you happen to use theoretical chemistry programs in your line of work then this blog will probably be of some interest to you and your colleges. During this development project I intend to create at least two scripts/programs which will accomplish the following tasks: Perform NEB (Nudged Elastic Band) simulations using Gaussian, and link this script to Jmol so as to create a seamless GUI that does not require a user to understand command line. Also, I intend for my code to be completely free and well supported. Lofty ambitions maybe, but we'll see what I can get done over the course of the summer.

If you do happen to have a vested interest in this project, please let me know! I am very open to feedback and I am excited to have input from any source!