#include <math.h>
#include "moldyn.h"
+#include "report/report.h"
int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
switch(type) {
case CUBIC:
- ret=cubic_init(a,b,c,lc,atom,NULL);
+ origin.x=0.5*lc;
+ origin.y=0.5*lc;
+ origin.z=0.5*lc;
+ ret=cubic_init(a,b,c,lc,atom,&origin);
break;
case FCC:
ret=fcc_init(a,b,c,lc,atom,NULL);
for(i=0;i<a;i++) {
r.y=o.y;
for(j=0;j<b;j++) {
+ r.z=o.z;
for(k=0;k<c;k++) {
- r.z=o.z;
v3_copy(&(atom[count].r),&r);
count+=1;
r.z+=lc;
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;
}
/* energy scaling factor */
energy_scale=moldyn->count*EV;
+printf("debug: %f\n",moldyn->atom[0].f.x);
/* calculate initial forces */
potential_force_calc(moldyn);
+printf("debug: %f\n",moldyn->atom[0].f.x);
/* some stupid checks before we actually start calculating bullshit */
if(moldyn->cutoff>0.5*moldyn->dim.x)
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,
if(!(i%v)) {
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);
+ printf("\rsched: %d, steps: %d, debug: %f",
+ sched->count,i,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 */