A Practical Introduction to the C Language for Computational Chemistry. Part 3

Sphere. From Space, from Space, Sir: whence else?

Square. Pardon me, my Lord, but is not your Lordship already in Space, your Lordship and his humble servant, even at this moment?

Sphere. Pooh! what do you know of Space? Define Space.

Square. Space, my Lord, is height and breadth indefinitely prolonged.

Sphere. Exactly: you see you do not even know what Space is. You think it is of Two Dimensions only; but I have come to announce to you a Third — height, breadth, and length.

Square. Your Lordship is pleased to be merry. We also speak of length and height, or breadth and thickness, thus denoting Two Dimensions by four names.

Sphere. But I mean not only three names, but Three Dimensions.

Adapted from: 
Flatland: A romance of many dimensions by Edwin A. Abbott


In the part 2 of this tutorial, we have learnt how to use arrays and how to read atomic coordinates from a file. In th eappendix, youc an find an example of solution of the exercises given in the previous tutorial.

In this third part, we are going to learn how to generate three-dimensional coordinated of atoms in a cubic crystal lattice and how to calculate non-bonded molecular potential and the force acting among them.

For this purpose, we need to modify the subroutines for the calculation of one dimnsional Lennard-Jones potentials and their forces developed in the previous tutorial to account to use directly the atomic coordinate of particles located in a tridimensional space.

The Cartesian distance of two points in a Euclidean space is given by:

r_{ij}=\sqrt{ (x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_i)^2}

where (x_i,y_i,z_i) and (x_j,y_j,z_j) are the coordinates of the two point particles. This expression will be used to modify the LJ energy potential and the components of the atomic forces.  The latter are derived taking in account the vectorial nature of the forces in the following components:

f(x_i)=(x_i-x_j)\frac{24\epsilon}{r_{ij}^2} \left[ 2\left(\frac{\sigma}{r_{ij}}\right)^{12} -  \left(\frac{\sigma}{r_{ij}}\right)^{6}  \right]

f(y_i)=(y_i-y_j)\frac{24\epsilon}{r_{ij}^2} \left[ 2\left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^{6} \right]

f(z_i)=(z_i-z_j)\frac{24\epsilon}{r_{ij}^2} \left[ 2\left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^{6} \right]

these expressions give the components of forces acting on the atom i. The force on the atom is given by \vec{f_j}=-\vec{f_i}, by the Netwon’s third law:

Lex. III. Actioni contrariam semper et equalem esse reactionem: sive corporum duorum actiones in se mutuo semper esse aequales et in partes contrarias dirigi.

Isaac Newton. Philosophiae Naturalis Principia Mathematica. London, 1686

Finally, the square root function in C is calculated by the sqrt().

EXERCISE 2: modify the program in the part 2 of this tutorial to include the possibility to read a file containing the coordinated in xyz format of two LJ particles. Using those coordinates, calculate the potential and forces acting between the two particles. Change the value of the coordinates and report the same value for different distances. 


Atoms in molecules are surrounded by an electron density that depend on by their tendency to delocalize their electron (these properties is usually referred to as electronegativity). The electron distribution produces and electric field around the molecule that can be matched by assigning partial electric point charges to each atom (I have explained in another article how this is accomplished). Therefore, we can add to the expression of the LJ potential a Coulomb potential term that describe the electrostatic pair interactions of charged atoms.

V(r_{ij})_{coul} = \frac{q_i q_j}{4\epsilon_0\epsilon_r r_{ij}}

The force components for the Coulomb term on the particle i interacting with the particle j are given by the expressions

f(x_i)_{coul}=\frac{q_i q_j}{4\pi\epsilon_0 \epsilon_r r_{ij}^2} \frac{(x_i-x_j)}{r_{ij}}

f(y_i)_{coul}=\frac{q_i q_j}{4\pi\epsilon_0\epsilon_r r_{ij}^2} \frac{(y_i-y_j)}{r_{ij}}

f(z_i)_{coul}=\frac{q_i q_j}{4\pi\epsilon_0\epsilon_r r_{ij}^2} \frac{(z_i-z_j)}{r_{ij}}

The q_i  and q_j are the partial charges and \epsilon_0 and \epsilon_r the dielectrioc constant of the vacuum and of the medium. The force on the atom is given by the Netwon’s third law. To extend the calculation to many-particles systems, we need to perform a double sum is used to calculate all the pair-wise interactions as follows

\begin{aligned}V(r_1,r_2,\cdots,r_n) &=  \\ &+ \sum_i \sum_{j>i}  4\epsilon_{ij}\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]\\  &+ \sum_i \sum_{j>i} \left(\frac{q_i q_j}{4\epsilon_0\epsilon_r r_{ij}} \right) \end{aligned}

The index i>j means that in the second sum we have to exclude all the pair already accounted by the first index. If we have four particles for example we need to calculate the following 4*(4-1)/2=6 pairs

i:   1    1    1   2    2     3    

j:   2     3   4   3    4     4

when the index i=2 the index j start from 3 since 2-1 has been calculated and 2-2 does not make sense.

To program this double summation, we need to use a two nested loops with two indexes (i and j). The following pseudo-code shows how to make the double loops. The potential and forces are updated using the unary operators += and they need to be initialized to zero. The results will be printed at the end of the calculation using a loop on the number of atoms (na).

Vt=0; /* total potential */
for (i=0;i<na;i++) { 
V[i] =0.0; } /* initialize the forces */
for (i=0;i<na;i++) {
   for (j=i+1; j< na; j++) {

calculate the distance using xx[i],xx[j],yy[i],yy[j],zz[i],zz[j]
calculate the LJ potential and forces

V+= V[i] + ….
Vt+= V[i] + ….
F[i]+= F[i] + …
/* Add here the printout of the potential for each atom and the total potential */

EXERCISE 3: implement the double loop to calculate the potential energy and forces for the coordinates generated by the program from EXERCISE 1.


Atomic coordinates for simple regular systems, such as a cubic crystal lattice, can be easily generated using nested for-loops. Let start with the example of a simple one-dimensional lattice, where the atoms are lined up along the x-axis. This means that coordinates of these atoms along the other two directions are null. The coordinate of a given atom along the x direction are generated by multiply the lattice index for the lattice size.  The following piece of code show how to implement the algorithm using for-loops. The variable xaxis is the lattice space in Angstrom units and nx is the number of lattice points along the x direction.

xaxis = 3.0; /* this value can be given in input */
nx = 5;
xx[0] = 0;
for (i=1;i<=nx;i++) {
  xx[i] = i*xaxis;

We can extend the method for a 2 D lattice using a nested for-loop as in this code snippet:

xaxis = 3.;
yaxis = 3.; 
xx[0] = 0.;
yy[0] = 0.;
for (i=0;i<nx;i++) {
   xx[i] = i*xaxis;
   for (j=0;j<ny;j++) {
       yy[j] = j*yaxis;

Now putting together what we have learnt from the last tutorial, we can write a code that generate the coordinates of a squared two-dimensional lattice and save them in xyz format. The code is reported in the appendix. The coordinated in xyz format (file name out.xyz) can be visualized them using the program rasmol. If you run from the comand line type rasmol –xyz  2dlattice.xyz.

With nx=10, and ny=1 you should get the the 1D-dimensional lattice that visualized with rasmol (background white, spacefill) will look like the following picture

With nx=10, and ny=10 you should get the the 2D-dimensional lattice that visualized with rasmol (background white, ball&stick) will look like the following picture

EXERCISE 1. Extend the 2D lattice program to calculate coordinates of atoms in a 3D lattice for a given number of lattice point along the 3 axis (nx, ny, nz) and for given lattice sizes (xaxis, yaxis, zaxis). You may want also to add the possibility to input these parameters from the keyboard to make the program more user-friendly.



Table I: Note the units express the final energies in kJ/mol and the distance in nanometers

PROGRAM I: Two dimensional cubic lattice generator

   PROGRAM NAME : 2D reticulum generator 
DESCRIPTION  : the program generate one or two dimensional  
lattices having nx,ny alternating atoms in the x 
and y direction, respectively. The atoms are separated
ati equal distances xsize,ysize
number of atoms in each dimension to be set
AUTHOR :  Danilo Roccatano 
DATE   : BREMEN, 2007

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <unistd.h> 

int  main() {
    char sym1[]="Na";
    char sym2[]="Cl";
    char csym[3];
    float xaxis,yaxis;
    int nx,ny;
    int nn,n,nsel;
    float xx[10];
    float yy[10];
    float zz;

     * Initialize variables


    /*Opening output file*/
    FILE *fout;
    if (!(fout=fopen("2dlattice.xyz","w"))){
        printf ("Error opening file %s\n", "2dlattice.xyz");
    /*Writing output header*/
    fprintf(fout,"%i %s%s-atoms in two dimensional rectangular lattice\n",n,sym1,sym2);
    /*Computing and writing the lattice coordinates*/
    for(int i=0;i<nx;i++){
        for(int j=0;j<ny;j++){
            if (nn%2) strcpy(csym, sym1);
            else strcpy(csym,sym2);
            fprintf(fout,"%2s %8.3f %8.3f %8.3f %5d\n", csym, xx[i], yy[j], zz,nn);
    printf("Generated a rectangular lattice of %i %s%s-atoms. \n",n,sym1,sym2);


This is an example of solution of the exercises for the last practical. The program read the parameters from the file “ljparam.dat” containing the following data.

LJ parameters
 Atom_name   sigma(A)       epsilon(K)
 Ar           3.336000      141.200000
 Kr           3.575000      191.400000
 Xe           3.924000      257.400000 

The program ask to choose two atoms types, then it calculates the combined parameters and ask for the maximum value of distance to calculate the value of the potential and the force. The values are stored in the file “Results.out“, and they can be plotted with xmagrace (see the graph at the end).

#include <stdio.h>
#include <math.h>
#include <string.h>
int main() {
    char line[160],title[160];
    char aty[10];
    char yn[2];
    int k=1;
    int ia1,ia2;
    int np; // number of parametes in the library
    float V,s,s1,s2,e,e1,e2,F;
    float x,rmin,xmin,xmax,xinc;
    float C6,C12;
    float s6,ix6,ix;
    float sig[10];
    float eps[10];
    float kb=8.3145e-3; /*Boltzmann constant*Avogadro Number  */
    FILE *fd;
    FILE *fout;


    if (!(fd=fopen("ljparam.dat","r"))){
        printf("\nError opening file %s\n", "ljparam.dat");

    fgets (title, 80, fd);
    fgets (line, 80, fd);
    sscanf (line,"%d",&np);
    fgets (line, 80, fd);
    printf (" %s  \n",line);
    while (fgets (line, 80, fd) != NULL){
        if (line[0] != '\n' && strlen(line) > 10) {
            sscanf (line,"%s%f%f",&aty[2*k],&sig[k],&eps[k]);
            printf ("%i %s  %8.3f   %8.3f\n",k, &aty[2*k],sig[k],eps[k]);

    if (!(fout=fopen("Results.out","w"))){
        printf ("Error opening file %s\n", "Results.out");

    /*Input for the atom types*/
    printf("Enter the atomtype for particle 1: ");
    printf("Enter the atomtype for particle 2: ");

    /*Berthelot-Lorenz mixing, etc.*/


    /*Output the results on the screen*/
    printf("Mixed Lennard-Jones parameters:\n");
    printf(" Epsilon:  %f kJ/mol \n",e);
    printf(" R min  :  %f nm     \n",rmin);
    printf(" Sigma  :  %f nm\n",s);
    printf(" C6     :  %f kJ*nm^6/mol \n",C6);
    printf(" C12    :  %f kJ*nm^12/mol \n\n",C12);
    /* Estimate the minimum distance to calculate the LJ function */
    printf(" Value of Potential V(xmin=%f nm):  %f kJ/mol \n",xmin,V);

    /*Output the results in the file  */
    fprintf(fout,"#Lennard-Jones parameters:\n\n");
    fprintf(fout,"# Atom Type 1  :  %c%c \n",aty[2*ia1],aty[2*ia1 + 1]);
    fprintf(fout,"#     Sigma    :  %f nm/10\n",sig[ia1]);
    fprintf(fout,"#   Epsilon    :  %f K\n\n",eps[ia1]);
    fprintf(fout,"# Atom Type 2  :  %c%c \n",aty[2*ia2],aty[2*ia2 + 1]);
    fprintf(fout,"#     Sigma    :  %f nm/10\n",sig[ia2]);
    fprintf(fout,"#   Epsilon    :  %f K\n\n",eps[ia2]);
    fprintf(fout,"# Mixed Sigma  :  %f nm\n",s);
    fprintf(fout,"# Mixed Epsilon:  %f kJ/mol \n",e);
    fprintf(fout,"#    C6        :  %f kJ*nm^6/mol \n",C6);
    fprintf(fout,"#   C12        :  %f kJ*nm^12/mol \n\n",C12);
    fprintf(fout,"# Distance [nm] Potential [kJ/mol] Force [nN]\n");

    /*Input for the value of the maximum distance range*/
    do {
        printf("\nEnter the value of the maximum distance distance (in nm): ");
        printf("\nEnter the distance increment (in nm): ");
    } while (strcmp(yn,"y"));

    for (x=xmin;x<=xmax;x+=xinc){ 
        /*Calculate the LJ potential and the force between xmin and xmax  x*/
        F=(12.*ix*C12*ix6*ix6-6.*ix*C6*ix6)/602.2; /*to obtain [nN]*/
        fprintf(fout," %f    %f    %f\n", x,V,F);
    return 0;

About Danilo Roccatano

I have a Doctorate in chemistry at the University of Roma “La Sapienza”. I led educational and research activities at different universities in Italy, The Netherlands, Germany and now in the UK. I am fascinated by the study of nature with theoretical models and computational. For years, my scientific research is focused on the study of molecular systems of biological interest using the technique of Molecular Dynamics simulation. I have developed a server (the link is in one of my post) for statistical analysis at the amino acid level of the effect of random mutations induced by random mutagenesis methods. I am also very active in the didactic activity in physical chemistry, computational chemistry, and molecular modeling. I have several other interests and hobbies as video/photography, robotics, computer vision, electronics, programming, microscopy, entomology, recreational mathematics and computational linguistics.
This entry was posted in Programming, Research, Science Topics. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.