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[count].mass);
92 v=sigma*rand_get_gauss(random);
94 p_total.x+=atom[count].mass*v;
96 v=sigma*rand_get_gauss(random);
98 p_total.x+=atom[count].mass*v;
100 v=sigma*rand_get_gauss(random);
102 p_total.x+=atom[count].mass*v;
105 /* zero total momentum */
106 for(i=0;i<count;i++) {
107 v3_scale(&delta,&p_total,1.0/atom[count].mass);
108 v3_sub(&(atom[count].v),&(atom[count].v),&delta);
111 /* velocity scaling */
112 scale_velocity(atom,count,t);
117 int scale_velocity(t_atom *atom,int count,double t) {
123 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
127 e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
128 c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
130 v3_scale(&(atom[count].v),&(atom[count].v),(1.0/c));
135 double get_e_kin(t_atom *atom,int count) {
143 e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));