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

Controlling complexity is the essence of computer programming.

Brian W. Kernighan. in Software Tools (with PJ Plauger), 1976.

THE FUNCTIONS

A C program is a collection of functions. A C function is equivalent to the subroutiing in FORTRAN or BASIC and procedures in PASCAL, PERL or PYTHON programming languages. It is a portion of the program that cannot be executed indipendently but only as part of another program. The function contains a specific algorithm or a stand alone procedure. You have already used several library functions in your previous programs. Output command for priting or reading files (such as printf(), openf()), mathematical functions (sqrt(), cos() are library or intrinsic functions as well. Other libraries functions, we can classify them as follows

  • Input/output functions. Input/output on the computer devices (e.g. output to the terminal, printer, hard disk, input from keyboard). It is usually used with #include <stdio.h>;
  • String manipulation functions. This library contains common operation on strings (e.g. concatenation, length, search and extraction of substrings). It is usually used with #include <string.h>;
  • Mathematical functions. Mathematical calculations (e.g. trigonometrics functions, exponentiation, square root extraction). It is usually used with #include <math.h>;
  • Graphical functions. Function for graphics operations (open a graphical window and canvas) and drawing graphical primitives (e.g. points, line, curves).
  • Operative system control functions. Operation requiriing allocation of the computer resources or devices (e.g. date and time, allocation of memory). It is usually used with, for example, #include <time.h>;
  • Data conversion functions. Operation for data conversion (e.g. change characters type, ascii to integer). It is usually used with #include <ctype.h>;

To use these function, you need to use the precompiler instruction #include at the beginning of the program. The compiler use by deafult the standard library #include <stdlib.h>;

You can write your own functions and it is called user defined functions. The use of function allow to structure the program and make easier its organization and reading. C language is structured around the use of functions. The function main is a function itself and it contains calls to other functions. both intrinsic and user defined functions.

The main function is the first function called by the operating system in your program. Every program must have exactly one main function. In principle, only code in main and functions called directly or indirectly from main will be executed. The main function may be implemented without arguments and has a return type of int:

int main () {

… // actual program code

 return 0;

}

The value returned from the main function is supplied to the operating system. As a standard, a value of 0 signals no error during program execution.

The definition of a function in C follows

  <return type>  function name ( <argument list> ) {

[ statments ] 

}

When a function does not return a result, then the return type is void. To return a value from a function, for this purpose the C language provides the keyword return

The value can be passed to the function by value or by reference. C passes parameters “by value” which means that the actual parameter values are copied into local storage. The caller and called functions do not shares any memory. This scheme is fine for many purposes, but it has two disadvantages.

  • Because the called function has its own copy, modifications to that memory are not communicated back to the caller. Therefore, value parameters do not allow the called function to communicate back to the caller. The function’s return value can communicate some information back to the caller, but not all problems can be solved with the single return value.
  • Sometimes it is undesirable to copy the value from the caller to the called function because the value is large and so copying it is expensive, or because at a conceptual level copying the value is undesirable.

Example

double square ( const double y ) {

  return y*y;

}  

The function is called using variable arguments and returns one or more values.

<function name> ( <argument1., <argument2>, . . . );

For example, one square is defined then you can pass the number 4.3 to this function by simply call square(4.3) and assing to the variable y the result:

double y = square( 4.3 );

The function is inserted before the definition of the main function or a a library file to be caalled using the #include precompiler statment. One can not define a function in the body of another function. A complete example of a program using the user defined function square is the following

#include <stdio.h>;

double square ( const double y ) { return y*y; }   

int main() {

    const double  m=10;

    double n=square(m);

    printf ("%d\n", n);

}

You can user function without other user defined functions. In the following example the function that calculate the cube of a number use the function square for its purpose.

#include <stdio.h>;

using namespace std;

double square ( const double y ) { return y*y; }   

double cube ( const double x ) {

    return square( x ) * x;

}

int main() {

    const double  m=10;

    double n=cube(m);

    printf ("%d\n", n);

}

The alternative is to pass the arguments “by reference”. Instead of passing a copy of a value from the caller to the called function, pass a pointer to the value. In this way there is only one copy of the value at any time, and the caller and called function both access that one value through pointers. Variables can be defined as global or local variables. Global variables can be accessed by all the function in the c program and are defined outside the function. Local variables are defined inside function. They are created when the function is used and then deleted from the memory as the function is exit.

As an example, we are going to modify the LJ Potential Calculator program shown in Part 2 to improve its structure by using functions.

We are going to move the code in the main(0 function in two function. One structure ReadParam() will be used to read the parameter from the file LJparam.dat and the other WriteOutput() to write the output on the Result.dat file. You can see that the modify main() function is semplified and more readable.

int main() {
    char aty[10];
    int ia1,ia2;
    int np; // number of parametes in the library
    float s,e;
    float C6,C12;
    float sig[10];
    float eps[10];
    float kb=8.3145e-3; /*Boltzmann constant*Avogadro Number  */
    float s6;
    float rmin;

    /* Read parameter file */

    ReadParam(np,aty,sig,eps);

    // Input for the atom types

    printf("\n");
    printf("Enter the atomtype for particle 1: ");
    scanf("%i",&ia1);
    printf("Enter the atomtype for particle 2: ");
    scanf("%i",&ia2);
    printf("\n");

    // Berthelot-Lorenz mixing

    s=0.5*(sig[ia1]+sig[ia2])/10.;
    e=kb*sqrt(eps[ia1]*eps[ia2]);
    s6=pow(s,6);
    C6=4*e*s6;
    C12=4*e*s6*s6;
    rmin=pow(2.,1./6.)*s;

    // 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);

    // Write the results on the file 

    WriteOutput(ia1,ia2,s,e,C12,C6,aty,sig,eps);

    return 0;
}

It is possible to futher semplify this part by moving in other functions the input of the parameters and the calculation and printing of the mixed parameters but for the moment shall we analyze these modifications. The LJ parameters for atom types are now read in the following ReadParam() function.

void ReadParam(int np,char *aty,float *sig,float *eps) {
    char line[160],title[160];
    int k;
    FILE *fd;

    /*
     * READ LJ PARAMETERS FROM A LIBRARY FILE
     */

    printf("\nREADING THE LJ PARAMETER FILE \n");
    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]);
            k++;
        }
    }
    fclose(fd);
}

The arrays with sigma (sig[]) and epsilon (eps[]) parameters as well as atom names (aty[]) are returned to the main program using pointers therefore, in the argument list of the function they are defined as pointers (using the *). In the function, the file is also opened and closed so the parameter of the file stream are assigned locally.

The choice of the atoms type to use for the two interacting atoms as well as the calculation and printing of the mixing parameters is left for the moment in the main() function. However, the saving on the file Result.dat of the output data and of the tabulated value of the function is moved in the function WriteOutput() function. Also in this case, pointers are used to pass the arrays to the function. The code in the function is the same as the original program with the exception the calculation of the xmin, the starting point of the tabulated function. We modify the program to calculate as starting minimum point the value of distance that give a value of the potential equal in magnitude to the value of the mixed parameter epsilon.

One significant numerical problem is the calculation of the roots of simple equations. Specifically, we want to find an algorithm capable of finding the numerical value of the unknown x that solves the algebraic problem F(x)=0.

In our case, we want to find a point x< sigma such that F(x)=VLJ(x)-epsilon=0. As shown in the Figure, by shifting the function VLJ(x) by epsilon, the point that we are looking for is a zero of the function. There are different numerical methods to obtain this result. As we can easily calculate the derivative of this function, a fast and easy-to-implement method is the Newton-Raphson one. The basic idea is to use the function’s derivative to calculate the function’s tangent at the starting point (indicated as x0 in Figure). Then the point of intersection of the tangent line with the X-axis (x1 in the Figure) is calculated. If the function is across the x-axis, the new point will give a first closer approximation of the zero for the function. Hence by calculating the tangent of the function in the new point again, a next point (x2) is obtained. As shown in the Figure, the new point is getting closer to the root of F(x)=0. That is because the absolute difference between |x_2-x_1| is less than |x_1-x_0|. By reiterating the process, it is possible to get closer and closer to the root and use the absolute difference between the current point with the previous one. It is also possible to define criteria to end the iteration when a given level of accuracy is reached. The implementation of the algorithm is relatively straightforward and it can be found in introduction textbooks of numerical analysis.

In the program, the Newton-Raphson was implemented in the function newton(). The functions f() and df() contains the shifted Lennard-Jones function and its derivative. The variable accuracy set the level of accuracy to calculate the root. In thi case it is set to 10^{-5). The newton() function is called in the WriteOutput() function before tabulation of the function to calculate the xmin value (the lower bound of the graph). The upper bound (xmax) and the x-increment are requested as input.

float  f(float e, float C12, float C6, double x) {
    float  ix6=1./pow(x,6.);
    float  V=C12*ix6*ix6-C6*ix6-e;
    return V;
}

float  df(float C12, float C6, double x) {
    float  ix=1./x;
    float  ix6=1./pow(x,6.);
    float  F=-(12.*ix*C12*ix6*ix6-6.*ix*C6*ix6);
    return F;
}

float  newton(float x1, float e, float C12,float C6) {

    /*
     * This function perform the search of zero of a function 
     * V(x)-eps=0 using the Newton method 
     */

    float   x, fx, fx1;
    float accuracy=1e-5;
    x1=0.02; // Start with a small value of x
    x = 0;
    while (fabs(x1 - x) >= accuracy)
        // Loop until the variation of x is less that the 
        // assigned accuracy 
    {
        x = x1;        // Assign the variable x1 equal to x
        fx = f(e,C12,C6,x);  // Calculate the V(x)-epsilon 
        fx1 = df(C12,C6,x); //Calculate value of f'(x) 
        x1 = x - (fx / fx1);    // Newton  iteration 
    };   

    return x1;
}

As last note, the following line of the codes at the beginningof the program

float  f(float, float , float , double);
float  df(float , float , double);
float  Newton(int, float);
void ReadParam(int ,char *,float *,float *);  
void WriteOutput(int ,int ,float ,float ,float ,float ,char *,float *,float *);

is the so-called function prototyping. The function prototype is a function declaration that specifies the function’s name and input/output interface of a function but omits the function body. Specifically, it provides the return type of the data that the function will return, the number, the order and the data type of arguments passed to the function. So it is used by the C compiler to check the correct call of functions in the program.

It is all for now, and in the following tutorial, we will see how the code can be further structured and how we can use the graphical library to make a plot for the curve directly on the screen.

Remember to express your interest in my tutorials by pressing the Like button, by sharing them or by adding a comment.

I wish you a very HAPPY NEW YEAR!

APPENDIX

The complete source code of the LJCalc program.

/* PROGRAM: LJCalc
 *
 * DESCRIPTION:
 * This is a simple LJ Potential Calculator.
 * The program read the LJ parameters of a 
 * list of atoms from the file ljparam.dat.
 * The program ask to select two atoms 
 * and it generate the mixing interaction
 * potential and a tabulated graph of it.
 * The results are saved in the file Results.out. 
 *
 *
 * VERSION: 1.1
 * AUTHOR: Danilo Roccatano * (c) 2017-2022
 */

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

/* 
 * Function prototyping
 */

float  f(float, float , float , double);
float  df(float , float , double);
float  Newton(int, float);
void ReadParam(int ,char *,float *,float *);  
void WriteOutput(int ,int ,float ,float ,float ,float ,char *,float *,float *);

/*
 * FUNCTIONS
 */

float  f(float e, float C12, float C6, double x) {
    float  ix6=1./pow(x,6.);
    float  V=C12*ix6*ix6-C6*ix6-e;
    return V;
}

float  df(float C12, float C6, double x) {
    float  ix=1./x;
    float  ix6=1./pow(x,6.);
    float  F=-(12.*ix*C12*ix6*ix6-6.*ix*C6*ix6);
    return F;
}

float  newton(float x1, float e, float C12,float C6) {

    /*
     * This function perform the search of zero of a function 
     * V(x)-eps=0 using the Newton method 
     */

    float   x, fx, fx1;
    float accuracy=1e-5;
    x1=0.02; // Start with a small value of x
    x = 0;
    while (fabs(x1 - x) >= accuracy)
        // Loop until the variation of x is less that the 
        // assigned accuracy 
    {
        x = x1;        // Assign the variable x1 equal to x
        fx = f(e,C12,C6,x);  // Calculate the V(x)-epsilon 
        fx1 = df(C12,C6,x); //Calculate value of f'(x) 
        x1 = x - (fx / fx1);    // Newton  iteration 
    };   

    return x1;
}

void ReadParam(int np,char *aty,float *sig,float *eps) {
    char line[160],title[160];
    int k;
    FILE *fd;

    /*
     * READ LJ PARAMETERS FROM A LIBRARY FILE
     */

    printf("\nREADING THE LJ PARAMETER FILE \n");
    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]);
            k++;
        }
    }
    fclose(fd);
}

void WriteOutput(int ia1, int ia2, float s, float e, float C12, float C6, char *aty, float *sig,float *eps) {

    char yn[2];
    float x,xmin,xmax,xinc;
    float s6,ix6,ix;
    float V,F;
    FILE *fout;

    /*
     *   OPEN THE OUTPUT FILE
     */

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

    /*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);

    /* Input for the value of the maximum distance range
     * For plotting the LJ function
     */

    do {
        printf("\nEnter the maximum distance distance (in nm): ");
        scanf("%f",&xmax);
        printf("\nEnter the distance increment (in nm): ");
        scanf("%f",&xinc);
        printf("Confirm?(y/n)\n");
        scanf("%s",yn);
    } while (strcmp(yn,"y"));


    /* Estimate the minimum distance to calculate the LJ function */
    xmin=newton(s,e,C12,C6);
    printf("Value of xmin that give V(xmin)=epsilon: %f nm \n",xmin);

    fprintf(fout,"# Distance [nm] Potential [kJ/mol] Force [nN]\n");
    for (x=xmin;x<=xmax;x+=xinc){ 
        /*Calculate the LJ potential and the force between xmin and xmax  x*/
        ix=1./x;
        ix6=1./pow(x,6.);
        V=C12*ix6*ix6-C6*ix6;
        F=(12.*ix*C12*ix6*ix6-6.*ix*C6*ix6)/602.2; /*to obtain [nN]*/
        fprintf(fout," %f    %f    %f\n", x,V,F);
    }
    fclose(fout);
}

/*
 * Main program
 *
 */

int main() {
    char aty[10];
    int ia1,ia2;
    int np; // number of parametes in the library
    float s,e;
    float C6,C12;
    float sig[10];
    float eps[10];
    float kb=8.3145e-3; /*Boltzmann constant*Avogadro Number  */
    float s6;
    float rmin;

    /* Read parameter file */

    ReadParam(np,aty,sig,eps);

    // Choice of atom types of the two interacting atoms.

    printf("\n");
    printf("Enter the atomtype for particle 1: ");
    scanf("%i",&ia1);
    printf("Enter the atomtype for particle 2: ");
    scanf("%i",&ia2);
    printf("\n");

    // Berthelot-Lorenz mixing

    s=0.5*(sig[ia1]+sig[ia2])/10.;
    e=kb*sqrt(eps[ia1]*eps[ia2]);
    s6=pow(s,6);
    C6=4*e*s6;
    C12=4*e*s6*s6;
    rmin=pow(2.,1./6.)*s;

    // 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);

    // Write the results on the file 

    WriteOutput(ia1,ia2,s,e,C12,C6,aty,sig,eps);

    return 0;
}

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 )

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.