From 4871747c5c848e5881bea7949a41ceb589263841 Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 24 Jan 2007 17:47:05 +0000 Subject: [PATCH] still pressure (+ fixes with cubic lattce and ho potential ...) --- moldyn.c | 33 +++++++++++++++------------------ moldyn.h | 2 -- report/report.h | 11 ++++------- sic.c | 18 +++++++++++------- 4 files changed, 30 insertions(+), 34 deletions(-) diff --git a/moldyn.c b/moldyn.c index 817a727..f564075 100644 --- a/moldyn.c +++ b/moldyn.c @@ -16,6 +16,7 @@ #include #include "moldyn.h" +#include "report/report.h" int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { @@ -590,18 +591,9 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) { double temperature_calc(t_moldyn *moldyn) { - double double_ekin; - int i; - t_atom *atom; - - atom=moldyn->atom; - - double_ekin=0; - for(i=0;icount;i++) - double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v)); + /* assume up to date kinetic energy, which is 3/2 N k_B T */ - /* kinetic energy = 3/2 N k_B T */ - moldyn->t=double_ekin/(3.0*K_BOLTZMANN*moldyn->count); + moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); return moldyn->t; } @@ -682,16 +674,21 @@ double pressure_calc(t_moldyn *moldyn) { double v; t_virial *virial; + /* + * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i ) + * + * virial = f_i r_i + */ + v=0.0; for(i=0;icount;i++) { virial=&(moldyn->atom[i].virial); v+=(virial->xx+virial->yy+virial->zz); } - v*=ONE_THIRD; -printf("kieck mal: %f %f %f\n",v,moldyn->count*K_BOLTZMANN*moldyn->t,v/moldyn->count); - moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+v; - moldyn->p/=moldyn->volume; + /* assume up to date kinetic energy */ + moldyn->p=2.0*moldyn->ekin+v; + moldyn->p/=(3.0*moldyn->volume); return moldyn->p; } @@ -1196,6 +1193,7 @@ int moldyn_integrate(t_moldyn *moldyn) { if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) scale_volume(moldyn); + update_e_kin(moldyn); temperature_calc(moldyn); pressure_calc(moldyn); //thermodynamic_pressure_calc(moldyn); @@ -1203,7 +1201,6 @@ int moldyn_integrate(t_moldyn *moldyn) { /* check for log & visualization */ if(e) { if(!(i%e)) - update_e_kin(moldyn); dprintf(moldyn->efd, "%f %f %f %f\n", moldyn->time,moldyn->ekin/energy_scale, @@ -1236,7 +1233,7 @@ int moldyn_integrate(t_moldyn *moldyn) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); printf("\rsched: %d, steps: %d, debug: %f | %f", - sched->count,i,moldyn->p/ATM,moldyn->debug/ATM); + sched->count,i,moldyn->p/ATM,moldyn->p/ATM); fflush(stdout); } } @@ -1580,7 +1577,7 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { d*=eps; v3_scale(&force,&distance,d); v3_add(&(aj->f),&(aj->f),&force); - v3_scale(&force,&distance,-1.0*d); /* f = - grad E */ + v3_scale(&force,&force,-1.0); /* f = - grad E */ v3_add(&(ai->f),&(ai->f),&force); virial_calc(ai,&force,&distance); virial_calc(aj,&force,&distance); /* f and d signe switched */ diff --git a/moldyn.h b/moldyn.h index 7329099..69cee96 100644 --- a/moldyn.h +++ b/moldyn.h @@ -12,8 +12,6 @@ #include "random/random.h" #include "list/list.h" -#include "report/report.h" - /* * * datatypes diff --git a/report/report.h b/report/report.h index 7fd3180..3868515 100644 --- a/report/report.h +++ b/report/report.h @@ -37,7 +37,7 @@ static char report_start[]="\ static char report_end[]="\ \\begin{figure}[!h]\n\ \\begin{center}\n\ -\\includegraphics[width=10cm]{energy.eps}\n\ +\\includegraphics[width=16cm]{energy.eps}\n\ \\caption{Kinetic, potential and total energy over time}\n\ \\end{center}\n\ \\end{figure}\n\ @@ -50,15 +50,12 @@ unset log \n\ unset label \n\ set xtic auto \n\ set ytic auto \n\ -set title 'Energy vs. time' \n\ +set title 'Energy per atom vs. time' \n\ set xlabel 'Time [fs]' \n\ -set ylabel 'Energy [eV]' \n\ -plot \"energy\" using 1:2 title 'Kinetic energy' with lines , \"energy\" using 1:3 title 'Potential energy' with lines , \"energy\" using 1:4 title 'Total energy' with lines \n\ -#set size 1.0, 0.6 \n\ +set ylabel 'Energy per atom [eV]' \n\ set terminal postscript eps enhanced color dashed lw 1 'Helvetica' 14 \n\ set output 'energy.eps' \n\ -replot\ +plot \"energy\" using 1:2 title 'Kinetic energy' with lines , \"energy\" using 1:3 title 'Potential energy' with lines , \"energy\" using 1:4 title 'Total energy' with lines\ "; - #endif diff --git a/sic.c b/sic.c index b5d4596..999d860 100644 --- a/sic.c +++ b/sic.c @@ -59,12 +59,13 @@ int main(int argc,char **argv) { //set_potential2b(&md,tersoff_mult_2bp,&tp); //set_potential2b_post(&md,tersoff_mult_post_2bp,&tp); //set_potential3b(&md,tersoff_mult_3bp,&tp); - set_potential2b(&md,lennard_jones,&lj); - //set_potential2b(&md,harmonic_oscillator,&ho); + //set_potential2b(&md,lennard_jones,&lj); + set_potential2b(&md,harmonic_oscillator,&ho); /* cutoff radius */ //set_cutoff(&md,TM_S_SI); - set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5)); + //set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5)); + set_cutoff(&md,1.1*LC_SI); /* * potential parameters @@ -78,7 +79,8 @@ int main(int argc,char **argv) { lj.uc=lj.epsilon4*(lj.sigma12/pow(md.cutoff,12.0)-lj.sigma6/pow(md.cutoff,6)); /* harmonic oscillator */ - ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; + //ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; + ho.equilibrium_distance=LC_SI; ho.spring_constant=.1; /* @@ -113,13 +115,14 @@ int main(int argc,char **argv) { tersoff_mult_complete_params(&tp); /* set (initial) dimensions of simulation volume */ - set_dim(&md,8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),TRUE); + //set_dim(&md,8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),TRUE); + set_dim(&md,8*LC_SI,8*LC_SI,8*LC_SI,TRUE); /* set periodic boundary conditions in all directions */ set_pbc(&md,TRUE,TRUE,TRUE); /* create the lattice / place atoms */ - create_lattice(&md,FCC,LC_SI*0.5*sqrt(1.5),SI,M_SI, + create_lattice(&md,CUBIC,LC_SI,SI,M_SI, //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, ATOM_ATTR_2BP|ATOM_ATTR_HB, @@ -144,7 +147,8 @@ int main(int argc,char **argv) { // &r,&v); /* setting a nearest neighbour distance for the moldyn checks */ - set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */ + //set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */ + set_nn_dist(&md,LC_SI); /* set temperature */ //set_temperature(&md,273.0+1410.0); -- 2.20.1