#include <math.h>
#include "moldyn.h"
+#include "report/report.h"
int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
double temperature_calc(t_moldyn *moldyn) {
- double double_ekin;
- int i;
- t_atom *atom;
-
- atom=moldyn->atom;
-
- double_ekin=0;
- for(i=0;i<moldyn->count;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;
}
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;i<moldyn->count;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;
}
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);
/* 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,
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);
}
}
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 */