return 0;
}
+int set_nn_dist(t_moldyn *moldyn,double dist) {
+
+ moldyn->nnd=dist;
+
+ return 0;
+}
+
int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
if(x)
}
moldyn->count=count;
+ printf("[moldyn] created lattice with %d atoms\n",count);
while(count) {
moldyn->atom[count-1].element=element;
}
/* velocity scaling */
- scale_velocity(moldyn);
+ scale_velocity(moldyn,VSCALE_INIT_EQUI);
return 0;
}
-int scale_velocity(t_moldyn *moldyn) {
+int scale_velocity(t_moldyn *moldyn,u8 type) {
int i;
- double e,c;
+ double e,scale;
t_atom *atom;
atom=moldyn->atom;
* - velocity scaling (E = 3/2 N k T), E: kinetic energy
*/
- if(moldyn->t==0.0) {
- printf("[moldyn] no velocity scaling for T = 0 K\n");
- return -1;
- }
-
e=0.0;
for(i=0;i<moldyn->count;i++)
e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
- c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
+ scale=(1.5*moldyn->count*K_BOLTZMANN*moldyn->t)/e;
+ if(type&VSCALE_INIT_EQUI) scale*=2.0; /* equipartition theorem */
+ scale=sqrt(scale);
for(i=0;i<moldyn->count;i++)
- v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
+ v3_scale(&(atom[i].v),&(atom[i].v),scale);
return 0;
}
t_3dvec p;
t_moldyn_schedule *schedule;
t_atom *atom;
-
int fd;
char fb[128];
+ double ds;
schedule=&(moldyn->schedule);
atom=moldyn->atom;
/* calculate initial forces */
potential_force_calc(moldyn);
+ /* do some checks before we actually start calculating bullshit */
+ if(moldyn->cutoff>0.5*moldyn->dim.x)
+ printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
+ if(moldyn->cutoff>0.5*moldyn->dim.y)
+ printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
+ if(moldyn->cutoff>0.5*moldyn->dim.z)
+ printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
+ ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
+ if(ds>0.05*moldyn->nnd)
+ printf("[moldyn] warning: forces too high / tau too small!\n");
+
/* zero absolute time */
moldyn->time=0.0;
moldyn->energy=0.0;
for(i=0;i<count;i++) {
-printf("BAR %d %d\n",i,count);
/* reset force */
v3_zero(&(atom[i].f));