#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
+#include <sys/time.h>
+#include <time.h>
#include <math.h>
#include "moldyn.h"
memset(moldyn,0,sizeof(t_moldyn));
+ moldyn->argc=argc;
+ moldyn->args=argv;
+
rand_init(&(moldyn->random),NULL,1);
moldyn->random.status|=RAND_STAT_VERBOSE;
atom[ret].brand=brand;
atom[ret].tag=count+ret;
check_per_bound(moldyn,&(atom[ret].r));
+ atom[ret].r_0=atom[ret].r;
}
/* update total system mass */
return ret;
}
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
+ t_3dvec *r,t_3dvec *v) {
+
+ t_atom *atom;
+ void *ptr;
+ int count;
+
+ atom=moldyn->atom;
+ count=(moldyn->count)++;
+
+ ptr=realloc(atom,(count+1)*sizeof(t_atom));
+ if(!ptr) {
+ perror("[moldyn] realloc (add atom)");
+ return -1;
+ }
+ moldyn->atom=ptr;
+
+ atom=moldyn->atom;
+ atom[count].r=*r;
+ atom[count].v=*v;
+ atom[count].element=element;
+ atom[count].mass=mass;
+ atom[count].brand=brand;
+ atom[count].tag=count;
+ atom[count].attr=attr;
+ check_per_bound(moldyn,&(atom[count].r));
+ atom[count].r_0=atom[count].r;
+
+ /* update total system mass */
+ total_mass_calc(moldyn);
+
+ return 0;
+}
+
/* cubic init */
int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
return count;
}
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
- t_3dvec *r,t_3dvec *v) {
-
- t_atom *atom;
- void *ptr;
- int count;
-
- atom=moldyn->atom;
- count=(moldyn->count)++;
-
- ptr=realloc(atom,(count+1)*sizeof(t_atom));
- if(!ptr) {
- perror("[moldyn] realloc (add atom)");
- return -1;
- }
- moldyn->atom=ptr;
-
- atom=moldyn->atom;
- atom[count].r=*r;
- atom[count].v=*v;
- atom[count].element=element;
- atom[count].mass=mass;
- atom[count].brand=brand;
- atom[count].tag=count;
- atom[count].attr=attr;
-
- /* update total system mass */
- total_mass_calc(moldyn);
-
- return 0;
-}
-
int destroy_atoms(t_moldyn *moldyn) {
if(moldyn->atom) free(moldyn->atom);
atom=moldyn->atom;
moldyn->ekin=0.0;
- for(i=0;i<moldyn->count;i++)
- moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+ for(i=0;i<moldyn->count;i++) {
+ atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+ moldyn->ekin+=atom[i].ekin;
+ }
return moldyn->ekin;
}
char dir[128];
double ds;
double energy_scale;
+ struct timeval t1,t2;
//double tp;
sched=&(moldyn->schedule);
moldyn->tau_square=moldyn->tau*moldyn->tau;
moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
- /* energy scaling factor */
- energy_scale=moldyn->count*EV;
+ /* get current time */
+ gettimeofday(&t1,NULL);
/* calculate initial forces */
potential_force_calc(moldyn);
moldyn->tau_square=moldyn->tau*moldyn->tau;
moldyn->time_steps=sched->runs[sched->count];
+ /* energy scaling factor (might change!) */
+ energy_scale=moldyn->count*EV;
+
/* integration according to schedule */
for(i=0;i<moldyn->time_steps;i++) {
/* check for log & visualization */
if(e) {
- if(!(i%e))
+ if(!(moldyn->total_steps%e))
dprintf(moldyn->efd,
"%f %f %f %f\n",
moldyn->time,moldyn->ekin/energy_scale,
get_total_energy(moldyn)/energy_scale);
}
if(m) {
- if(!(i%m)) {
+ if(!(moldyn->total_steps%m)) {
momentum=get_total_p(moldyn);
dprintf(moldyn->mfd,
"%f %f %f %f %f\n",moldyn->time,
}
}
if(p) {
- if(!(i%p)) {
+ if(!(moldyn->total_steps%p)) {
dprintf(moldyn->pfd,
"%f %f %f %f %f\n",moldyn->time,
moldyn->p/BAR,moldyn->p_avg/BAR,
}
}
if(t) {
- if(!(i%t)) {
+ if(!(moldyn->total_steps%t)) {
dprintf(moldyn->tfd,
"%f %f %f\n",
moldyn->time,moldyn->t,moldyn->t_avg);
}
}
if(s) {
- if(!(i%s)) {
+ if(!(moldyn->total_steps%s)) {
snprintf(dir,128,"%s/s-%07.f.save",
moldyn->vlsdir,moldyn->time);
- fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT);
+ fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
+ S_IRUSR|S_IWUSR);
if(fd<0) perror("[moldyn] save fd open");
else {
write(fd,moldyn,sizeof(t_moldyn));
}
}
if(v) {
- if(!(i%v)) {
+ if(!(moldyn->total_steps%v)) {
visual_atoms(&(moldyn->vis),moldyn->time,
moldyn->atom,moldyn->count);
}
}
/* display progress */
- if(!(i%10)) {
- printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f",
+ if(!(moldyn->total_steps%10)) {
+ /* get current time */
+ gettimeofday(&t2,NULL);
+
+ printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
sched->count,i,
moldyn->t,moldyn->t_avg,
- moldyn->p_avg/BAR,moldyn->p/BAR,
- moldyn->volume);
+ moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
+ moldyn->volume,
+ (int)(t2.tv_sec-t1.tv_sec));
fflush(stdout);
+
+ /* copy over time */
+ t1=t2;
}
/* increase absolute time */
return 0;
}
+/*
+ * restore function
+ */
+
+int moldyn_load(t_moldyn *moldyn) {
+
+ // later ...
+
+ return 0;
+}
+
/*
* post processing functions
*/
}
}
+int analyze_bonds(t_moldyn *moldyn) {
+
+
+
+
+ return 0;
+}
+