2 * moldyn.c - molecular dynamics library main file
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
14 #include "math/math.h"
15 #include "init/init.h"
16 #include "random/random.h"
19 int create_lattice(unsigned char type,int element,double mass,double lc,
20 int a,int b,int c,t_atom **atom) {
28 if(type==FCC) count*=4;
29 if(type==DIAMOND) count*=8;
31 *atom=malloc(count*sizeof(t_atom));
33 perror("malloc (atoms)");
41 ret=fcc_init(a,b,c,lc,*atom,&origin);
44 ret=diamond_init(a,b,c,lc,*atom,&origin);
47 printf("unknown lattice type (%02x)\n",type);
53 printf("ok, there is something wrong ...\n");
54 printf("calculated -> %d atoms\n",count);
55 printf("created -> %d atoms\n",ret);
60 (*atom)[count-1].element=element;
61 (*atom)[count-1].mass=mass;
68 int destroy_lattice(t_atom *atom) {
75 int thermal_init(t_atom *atom,t_random *random,int count,double t) {
78 * - gaussian distribution of velocities
79 * - zero total momentum
80 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
85 t_3dvec p_total,delta;
87 /* gaussian distribution of velocities */
89 for(i=0;i<count;i++) {
90 sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
92 v=sigma*rand_get_gauss(random);
94 p_total.x+=atom[i].mass*v;
96 v=sigma*rand_get_gauss(random);
98 p_total.y+=atom[i].mass*v;
100 v=sigma*rand_get_gauss(random);
102 p_total.z+=atom[i].mass*v;
105 /* zero total momentum */
106 v3_scale(&p_total,&p_total,1.0/count);
107 for(i=0;i<count;i++) {
108 v3_scale(&delta,&p_total,1.0/atom[i].mass);
109 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
112 /* velocity scaling */
113 scale_velocity(atom,count,t);
118 int scale_velocity(t_atom *atom,int count,double t) {
124 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
128 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
129 c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
131 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
136 double get_e_kin(t_atom *atom,int count) {
143 for(i=0;i<count;i++) {
144 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
150 double get_e_pot(t_moldyn *moldyn) {
152 return(moldyn->potential(moldyn));
155 double get_total_energy(t_moldyn *moldyn) {
159 e=get_e_kin(moldyn->atom,moldyn->count);
160 e+=get_e_pot(moldyn);
165 t_3dvec get_total_p(t_atom *atom, int count) {
171 for(i=0;i<count;i++) {
172 v3_scale(&p,&(atom[i].v),atom[i].mass);
173 v3_add(&p_total,&p_total,&p);
182 * potentials & corresponding forces
186 /* lennard jones potential & force for one sort of atoms */
188 double potential_lennard_jones(t_moldyn *moldyn) {
197 double eps,sig6,sig12;
199 params=moldyn->pot_params;
204 sig12=params->sigma12;
207 for(i=0;i<count;i++) {
209 v3_sub(&distance,&(atom[j].r),&(atom[i].r));
210 d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
211 help=d*d; /* 1/r^4 */
213 d=help*help; /* 1/r^12 */
214 u+=eps*(sig12*d-sig6*help);
221 int force_lennard_jones(t_moldyn *moldyn) {
229 double eps,sig6,sig12;
233 params=moldyn->pot_params;
236 sig12=params->sigma12;
238 for(i=0;i<count;i++) v3_zero(&(atom[i].f));
240 for(i=0;i<count;i++) {
242 v3_sub(&distance,&(atom[j].r),&(atom[i].r));
243 v3_per_bound(&distance,&(moldyn->dim));
244 d=v3_absolute_square(&distance);
245 if(d<=moldyn->cutoff_square) {
246 h1=1.0/d; /* 1/r^2 */
255 v3_scale(&force,&distance,d);
256 v3_add(&(atom[j].f),&(atom[j].f),&force);
257 v3_sub(&(atom[i].f),&(atom[i].f),&force);