+double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
+
+ double tau;
+
+ tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
+ tau*=1.0E-9;
+ if(tau<moldyn->tau)
+ printf("[moldyn] warning: time step (%f > %.15f)\n",
+ moldyn->tau,tau);
+
+ return tau;
+}
+
+
+/*
+ *
+ * 'integration of newtons equation' - algorithms
+ *
+ */
+
+/* start the integration */
+
+int moldyn_integrate(t_moldyn *moldyn) {
+
+ int i;
+ int write;
+
+ write=moldyn->write;
+
+ /* calculate initial forces */
+ moldyn->force(moldyn);
+
+ for(i=0;i<moldyn->time_steps;i++) {
+ /* integration step */
+ moldyn->integrate(moldyn);
+
+ /* check for visualiziation */
+ if(!(i%write)) {
+ visual_atoms(moldyn->visual,i*moldyn->tau,
+ moldyn->atom,moldyn->count);
+ }
+ }
+
+ return 0;
+}
+
+/* velocity verlet */
+
+int velocity_verlet(t_moldyn *moldyn) {
+
+ int i,count;
+ double tau,tau_square;
+ t_3dvec delta;
+ t_atom *atom;
+
+ atom=moldyn->atom;
+ count=moldyn->count;
+ tau=moldyn->tau;
+
+ tau_square=tau*tau;
+
+ for(i=0;i<count;i++) {
+ /* new positions */
+ v3_scale(&delta,&(atom[i].v),tau);
+ v3_add(&(atom[i].r),&(atom[i].r),&delta);
+ v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
+ v3_add(&(atom[i].r),&(atom[i].r),&delta);
+ v3_per_bound(&(atom[i].r),&(moldyn->dim));
+
+ /* velocities */
+ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
+ v3_add(&(atom[i].v),&(atom[i].v),&delta);
+ }
+
+ /* forces depending on chosen potential */
+ moldyn->force(moldyn);
+
+ for(i=0;i<count;i++) {
+ /* again velocities */
+ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
+ v3_add(&(atom[i].v),&(atom[i].v),&delta);
+ }
+
+ return 0;
+}
+