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 subroutine in FORTRAN or BASIC and procedures in PASCAL, PERL, or PYTHON programming languages. A portion of the program cannot be executed independently 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 commands for printing or reading files (such as printf(), openf()), mathematical functions (sqrt(), cos() are a library or intrinsic functions as well. Other library functions, we can classify 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 operations 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 functions, and it is called user-defined functions. The use of function allows the program to structure and makes its organization and reading easier. C language is structured around the use of functions. The main() function is a function that 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,” meaning the actual parameter values are copied into local storage. The caller and called functions do not shares any memory. This scheme is acceptable for many purposes, but it has two disadvantages.
- Because the called function has its 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 to the caller, but not all problems can be solved with a single return value.
- Sometimes it is undesirable to copy the value from the caller to the called function because the value is considerable 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 library file to be called using the #include precompiler statement. 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 use functions without other user-defined functions. In the following example, the function that calculates the cube of a number uses 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. All the functions in the c program can access global variables and are defined outside the function. Local variables are defined inside the function. They are created when the function is used and then deleted from the memory as it exits.
For example, we will modify the LJ Potential Calculator program shown in Part 2 to improve its structure using functions.
We will move the code in the main() function into two functions. 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 modified main() function is simplified 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 further simplify this part by moving the input of the parameters into other functions and the calculation and printing of the mixed parameters. Still, for now, we shall 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 is assigned locally.
The choice of the atom types to use for the two interacting atoms and 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 the tabulated value of the function is now performed in the 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 except for 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 gives 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 
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 
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 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;
}
daniloroccatano.blog 