The merit of using ab-initio calculations rests in the ability to predict phenomenology associated with chemical processes, energies, structures, thermodynamic properties, and other properties of atomic and molecular systems that is difficult to determine through experimental means. In this thesis, I address the structure and growth of clusters of calcium and calcium carbonate molecules by using density functional hybrid methods within the framework of the general gradient approximation, B3PW91/6-311g(d). Additionally I extended the study to the reaction paths and activation energy in the degradation of GB by H2O and OH.
In Chapter 3 I have presented an analysis of small
calcium clusters containing up to 13 calcium atoms using the method/basis set
B3PW91/6-311g(d). I have determined the
structure, binding energies, and harmonic frequencies of vibration for several
isomers at each cluster size up to Ca13. The lowest-energy structures of Ca6 and Ca10
through Ca13 had not been reported before as global geometries for
calcium clusters. I predict that there
are two paths of cluster growth. At almost zero Kelvin, clusters grow by
twinning pentagonal bipyramids as the building block, while at temperatures closely
above room temperature, clusters grow building shells of atoms around a central
atom with local cubic symmetry. For Ca12, instead of the
icosohedron-minus-one-atom expected structure, we found that this structure is
not reachable within the temperature range studied here (below 500 K). The Ca12 lowest-energy isomer is
a structure with two twinned pentagonal bipyramids that has a decorated
face. I also found that this isomer
undergoes a structural transition at about 317 K and becomes an isomer with
local cubic symmetry elements as depicted in Fig. 3.15
(Cs because of slight distortions). The metallic character of calcium is far
from being attained at Ca13 because the HOMO-LUMO energy difference
is 0.8 eV, a value that is very large when compared to the small gap of 0.03 eV
observed in bulk fcc calcium at the edge of the filled band.
In Chapter 4, I investigated the reaction of GB with the
water molecule and the OH radical and I predicted several reaction paths,
transition structures, and activation energies. I found that the reaction of GB with OH with the formation of an
intermediate structure is a very difficult and an improbable reaction under
normal conditions because the initial barrier of 76.8 kcal/mole is needed to
start the reaction. It is observed that
if 76.8 kcal/mole are available, the formation of the (GB+OH)* neutral radical
is possible. If the neutral radical is
formed then two possible unimolecular reactions can proceed as shown in Fig.
4.14 and 4.15. The energy barrier that
is needed for the fluorine abstraction reaction is 25.6 kcal/mole, whereas only
20.8 kcal/mole are required to liberate HF.
In the reaction of GB and H2O, the
energy needed to overcome the barrier is 38.6 kcal/mol if calculated within the
GGA and 60.6 in the HF/D95** approximation.
In this type of reaction the products HF and IMPA are produced, having a
net relative energy of 6.69 kcal/mole.
The half-life of a solution of GB and water is approximately 100
hours. If a salt is added to the water
to increase the pH level (becoming more basic) the rate of degradation of GB
will increase. A solution of sodium
hydroxide (NaOH) is commonly used in the neutralization or destruction of GB
and other nerve agents.
In the study of clusters of CaCO3, I used the results of the density functional hybrid method to make adjustments in the parameters of the RIM potential. I located low minimum energy structure containing up to 500 atoms (100 molecules of CaCO3), studied the relationship between the binding energy and cluster size, and investigated the rate of diffusion as a function of temperature for different cluster sizes.
In the study of molecular processes and ab-initio calculations, there is room for additional future research. Future work include:
1. Find a semi-empirical model that describes the interactions between calcium atoms in the cluster based on ab-initio results and extend the study to larger structures.
2. Calculate the reaction rates associated with the transition structures and reactions that have been found in my study of the degradation of GB.
3. Study the reactions and transition structures of GB with basic or acidic compounds. This would lead to an understanding of how the degradation of GB is dependant on the ph-level.
4. Investigate the formation processes and the dynamics that lead to the formation of calcite and aragonite. The crystalline structure of calcite and aragonite contain millions of atoms and is formed over long periods of time. To study such processes one would need to perform large-scale simulations on a large number of atoms over a long time period. Such a study would require a larger and more powerful computer system.
The application software that was used in the research included Gaussian98W (for PC) and Gaussian98 (for UNIX platforms), Microsoft Visual C++, Borland C++, Netscape 4.7X, Rasmol, Chimes, and Paintshop Pro. Gaussian98 was used to perform the ab-initio energy calculations, geometric optimizations, population analysis, frequency calculations, reaction kinetics, optimization of transition structures, and Intrinsic Reaction Calculation (IRC). The information from the generated log files was used in much of the analysis. The plots were generated using Matlab routines.
The visualizations of the atom/molecular structures were generated by writing a C++ program that would take the coordinate information from the Gaussian98 log file, and save the coordinates in either ‘pdb’ or ‘xyz’ formats. Rasmol or Netscape, with the Chimes plug-in, could then be used in visualizing the 3D structure.
|
// Program for generating ‘pdb’ files from
Gaussian98 output files void
pdb1(void) { int i; FILE *fp; if
((fp=fopen("ts1.pdb","w"))==NULL) { printf("Cannot open
file\n"); exit(1); } fprintf(fp,"COMPND Calcium\n"); for (i=0;i<n_atoms;i++) { fprintf(fp,"ATOM "); fprintf(fp,"%d",i); if (i<10) fprintf(fp," "); if (i>9) fprintf(fp," "); if(element[i]==1) fprintf(fp,"H 1 "); if(element[i]==6) fprintf(fp,"C 1 "); if(element[i]==8) fprintf(fp,"O 1 "); if(element[i]==9) fprintf(fp,"F 1 "); if(element[i]==15)
fprintf(fp,"P 1 "); if(element[i]==16)
fprintf(fp,"S 1 "); if(element[i]==7) fprintf(fp,"N 1 "); if(element[i]==11)
fprintf(fp,"Na 1 "); if(element[i]==17)
fprintf(fp,"Cl 1 "); if(element[i]==20)
fprintf(fp,"Ca 1 "); fprintf(fp,"%6.3lf %6.3lf
%6.3lf\n",x[i],y[i],z[i]); } fclose(fp); } |
The ‘element’ array stores the atomic number of the atom and the arrays ‘x’, ‘y’, and ‘z’ contain the coordinates of each atom.
Animation sequences were accomplished by generating a series of frames in the ‘xyz’ format. A C++ program was written that would create an animation sequence in ‘xyz’ format from a molecular dynamics simulation. The animation of the vibrational modes was accomplished by using the coordinates from the Gaussian98 log file, the amplitudes of the coefficients of the eigenvectors, creating a sinusoidal oscillation with a frequency corresponding to the frequency of that vibrational normal mode. The ‘element’ array is the atomic number of the atom, the arrays ‘x’, ‘y’, and ‘z’ contain the coordinates of each atom, and the arrays ‘dx’, ‘dy’, ‘dz’ contain the amplitudes of the oscillation of the vibrational mode.
|
if ((fp=fopen(file,"w"))==NULL)
{ printf("Cannot open
file\n"); exit(1); } ic=0; nfreq=int(30.0); for(j=0;j<nfreq;j++) { ic=ic+1; t=2.0*(double)(j)*pi/nfreq; ct=cos(t); if (ic==1)
fprintf(fp,"%d\n-\n",n_atoms); if (ic==2)
fprintf(fp,"%d\n\\\n",n_atoms); if (ic==3)
fprintf(fp,"%d\n|\n",n_atoms); if (ic==4) { fprintf(fp,"%d\n\/\n",n_atoms); ic=0; } for(k=0;k<n_atoms;k++) { if(element[k]==1)
fprintf(fp,"H"); if(element[k]==6)
fprintf(fp,"C"); if(element[k]==8)
fprintf(fp,"O"); if(element[k]==9)
fprintf(fp,"F"); if(element[k]==15) fprintf(fp,"P"); fprintf(fp," %10.5lf %10.5lf
%10.5lf\n", (x[k]+dx[k]*ct),(y[k]+dy[k]*ct), (z[k]+dz[k]*ct)); }
} fclose(fp); } } |
In visualizing the animation with Netscape and the Chimes plug-in, a small chimes/rasmol script was written within the “embed” command in HTML.
|
<!—This
is an example of the script--> <embed
align=abscenter src="calcite/animation/anim30.xyz" bgcolor=#fcfcfc
frank=false script='zoom 140; wireframe 20; cpk 80; rotate z 0; rotate y 0' animfps=20
startanim=true width=350 height=350 options3d=specular animmode="loop"> |
During the course of the study several computer programs were used, some were routines I wrote and others were routines from books. Here is a list of the computer programs used in the thesis:
1. A good method for finding lowest energy reaction paths is the Nudged-Elastic Band method. This was used in conjunction with Gaussian98 to locate minimum energy paths of the reaction with GB and OH. I wrote a C++ program that implemented the NEB algorithm. The program would read the atomic coordinates and the force values that were generated from the Gaussian98 log file and use these value in locating the minimum energy reaction path. See section 2.6 for further information. Here is the routine that performs the NEB method
|
void
neb(int images)
(1) { double xx[30],yy[30],zz[30]; double dx[30],dy[30],dz[30]; double
fsx,fsy,fsz,true_forcex,true_forcey,true_forcez; double
x1,x2,y1,y2,z1,z2,tx,ty,tz,d1,d2,d3;
double dx1,dy1,dz1,dx2,dy2,dz2,cosi,fsw; double c1,c2,tn,tnn,kk; int i; // for (i=0;i<n_atoms;i++) { xx[i]=x[images][i]; yy[i]=y[images][i]; zz[i]=z[images][i]; } // // Forces in the x, y, and z direction // for(i=0;i<n_atoms;i++) { dx[i]=fx[images-1][i]; dy[i]=fy[images-1][i]; dz[i]=fz[images-1][i]; } // // Defines Switching function // dx1=0.0; dy1=0.0; dz1=0.0; dx2=0.0; dy2=0.0; dz2=0.0; d1=0.0; d2=0.0; d3=0.0; for
(i=0;i<n_atoms;i++) { x1=xx[i]-x[images-1][i]; y1=yy[i]-y[images-1][i]; z1=zz[i]-z[images-1][i]; x2=x[images+1][i]-xx[i]; y2=y[images+1][i]-yy[i]; z2=z[images+1][i]-zz[i]; dx1=dx1+x1; dy1=dy1+y1; // (2) dz1=dz1+z1; dx2=dx2+x2; dy2=dy2+y2; dz2=dz2+z2; d1=d1+x1*x1+y1*y1+z1*z1; d2=d2+x2*x2+y2*y2+z2*z2; } cosi=(dx1*dx2+dy1*dy2+dz1*dz2)/sqrt(d1*d2); fsw=0.5*(1.0+cos(pi*cosi)); // //
Calculates a vector that is tangent to the image // d1=sqrt(d1); d2=sqrt(d2); d3=sqrt(d3); tn=0.0; c1=0.0; c2=0.0; kk=1.0; for(i=0;i<n_atoms;i++) { x1=xx[i]-x[images-1][i]; y1=yy[i]-y[images-1][i]; z1=zz[i]-z[images-1][i]; x2=x[images+1][i]-xx[i]; y2=y[images+1][i]-yy[i]; z2=z[images+1][i]-zz[i]; tx=x1/d1+x2/d2; ty=y1/d1+y2/d2; tz=z1/d1+z2/d2; tn=tn+tx*tx+ty*ty+tz*tz; } |
tnn=sqrt(tn); (3) for(i=0;i<n_atoms;i++) { x1=xx[i]-x[images-1][i]; y1=yy[i]-y[images-1][i]; z1=zz[i]-z[images-1][i]; x2=x[images+1][i]-xx[i]; y2=y[images+1][i]-yy[i]; z2=z[images+1][i]-zz[i]; fsx=kk*(x2-x1); fsy=kk*(y2-y1); fsz=kk*(z2-z1); tx=(x1/d1+x2/d2)/tnn; ty=(y1/d1+y2/d2)/tnn; tz=(z1/d1+z2/d2)/tnn; c2=c2+fsx*tx+fsy*ty+fsz*tz; c1=c1+dx[i]*tx+dy[i]*ty+dz[i]*tz; } // // Begin calculation of force for each atom // for(i=0;i<n_atoms;i++) { x1=xx[i]-x[images-1][i]; y1=yy[i]-y[images-1][i]; z1=zz[i]-z[images-1][i];// x2=x[images+1][i]-xx[i]; y2=y[images+1][i]-yy[i]; z2=z[images+1][i]-zz[i]; tx=(x1/d1+x2/d2)/tnn; ty=(y1/d1+y2/d2)/tnn; tz=(z1/d1+z2/d2)/tnn; // // Calculates normal component of the true
force // true_forcex=dx[i]-c1*tx; true_forcey=dy[i]-c1*ty; true_forcez=dz[i]-c1*tz; // // Calculates tangential component (4) //of
the spring force // force = dot product of rtf . rt using the
images image-1, image, image+1 // fsx=kk*(dx2-dx1); fsy=kk*(dy2-dy1); fsz=kk*(dz2-dz1); // //
Defines total force on each atom of the image // fxx[i]=(true_forcex+c2*tx+fsw*(fsx-c2*tx))*1.6e-19*1.0e10;
fyy[i]=(true_forcey+c2*ty+fsw*(fsy-c2*ty))*1.6e-19*1.0e10;
fzz[i]=(true_forcez+c2*tz+fsw*(fsz-c2*tz))*1.6e-19*1.0e10; } } |
The coordinates of every atom of each image are located in the 3D array ‘x’, ‘y’, and ‘z’ and the forces acting on every atom of each image is stored in ‘fx’, ‘fy’, and ‘fz’. In the routine the ‘NEB’ forces being calculated are the forces acting on each atom in one ‘image’. The forces and the coordinates acting on every atom of each ‘image’ are the input to the routine. Once the ‘NEB’ force is calculated on the image, the atoms in the image are moved in a direction and magnitude depending on the force calculated in the ‘NEB’ routine. The variable ’j’ is the image number and the variable ‘i’ is the atom number.
|
for (j=1;j<np;j++) { neb(j); // // Position and velocity // for (i=0;i<n_atoms;i++) { if (element[i]==1) mass=mass_h; if (element[i]==6) mass=mass_c; if (element[i]==8) mass=mass_o; if (element[i]==9) mass=mass_f; if (element[i]==15) mass=mass_p; x[j][i] = x[j][i] +
fxx[i]/2.0/mass*dt*dt; y[j][i] = y[j][i] +
fyy[i]/2.0/mass*dt*dt; z[j][i] = z[j][i] +
fzz[i]/2.0/mass*dt*dt; }
} |
2. During the study of calcium and CaCO3 clusters I made use of the Fletcher-Reeves-Polak-Ribierre minimization as described in Numerical Recipes [67] to locate minimum energy structures using the MM potential for calcium clusters and the RIM potential for CaCO3 clusters.
3. In finding the best fit of the parameters from the RIM potential to the Gaussian98 calculations, I used a non-linear least square fitting employing the “crude” method and Newton’s method as described in Devries [69]. Refer to section 2.7.
4. Another program was written in C++ that performed translation and rotational about a coordinate system on GB to align the P-C bond along the z-axis and orient the P-F bond with the x-axis. This allowed for a better visualization and manipulation of the atoms in locating the transition structure of GB and OH.
5. In the study of calcium clusters a C++ program was written to calculate the thermodynamic properties of calcium clusters in the harmonic approximation as described in section 3.4.4. The program was useful in locating a structural transition in Ca12. The transition between structures occurred at approximately 317 Kelvin.
6. In the study of the CaCO3, I had to locate minimum energy structure for structures containing up to 500 atoms, or 100 CaCO3 molecules. A C++ program was written to perform the simulated annealing with Molecular Dynamics, as described in section 5.7. In the program I made use of the Velocity-Verlet algorithm. The input is the initial atomic coordinates and the initial velocity of the atoms. The atomic position and velocities are change with each time step (10-15 second) according to the atomic forces acting on each atom, which are derived from the RIM potential.
|
force(rx0,ry0,rz0); for
(it=1;it<ntime;it++) { printf("%d\n",it); icount=icount+1; ic=ic+1; icc=icc+1; // // Position and velocity // for
(i=0;i<n_atoms;i++) { if (element[i]==20) value=dt/2.0/mass_ca; if (element[i]==6) value=dt/2.0/mass_c; if (element[i]==8) value=dt/2.0/mass_o; rx[i] = rx[i] + vx[i]*dt + fx[i]*value*dt; ry[i] = ry[i] + vy[i]*dt + fy[i]*value*dt; rz[i] = rz[i] + vz[i]*dt + fz[i]*value*dt; vx[i] = vx[i] + fx[i]*value; vy[i] = vy[i] + fy[i]*value; vz[i] = vz[i] + fz[i]*value; rxt[i]=rx[i]*1.0e10; ryt[i]=ry[i]*1.0e10; rzt[i]=rz[i]*1.0e10; } // //
Calculate the Force // force(rxt,ryt,rzt); // //
Velocity // for
(i=0;i<n_atoms;i++) { if (element[i]==20)
value=dt/2.0/mass_ca; if (element[i]==6) value=dt/2.0/mass_c; if (element[i]==8) value=dt/2.0/mass_o; vx[i] = vx[i] + fx[i]*value; vy[i] = vy[i] + fy[i]*value; vz[i] = vz[i] + fz[i]*value; } |
// // Scaling // if
(ic==ncut) { ke=0.0; for
(i=0;i<n_atoms;i++) { if
(element[i]==20) mass=mass_ca; if
(element[i]==6) mass=mass_c; if
(element[i]==8) mass=mass_o; temp[3*i]=rxt[i]; temp[3*i+1]=ryt[i]; temp[3*i+2]=rzt[i]; ke=ke+0.5*mass*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]); } ke1=ke/charge; rim1=rim(temp); e1=rim1+ke1; temperature=ke1*2.0*11604.45/3.0/n_atoms; fprintf(fp,"%d %lf
%lf %lf %lf\n", it,rim1,ke1,e1,temperature); printf("%d %lf
%lf %lf %lf\n",
it,rim1,ke1,e1,temperature); if (temperature<400) scale=1.01; else scale=1.0; if
(temperature>600) scale=0.80; for
(i=0;i<n_atoms;i++) { vx[i]=vx[i]*scale; vy[i]=vy[i]*scale; vz[i]=vz[i]*scale; } ic=0; }
if (icount==icut) { icount=0; itt=itt+1; anim(rxt,ryt,rzt,itt); if (itt>2) itt=-1; }
} |