return 0;
}
+int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) {
+
+ moldyn->bondlen[0]=b0*b0;
+ moldyn->bondlen[1]=b1*b1;
+ if(bm<0)
+ moldyn->bondlen[2]=b0*b1;
+ else
+ moldyn->bondlen[2]=bm*bm;
+
+ return 0;
+}
+
int set_temperature(t_moldyn *moldyn,double t_ref) {
moldyn->t_ref=t_ref;
break;
case VISUAL_STEP:
moldyn->vwrite=timer;
- ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
+ ret=visual_init(moldyn,moldyn->vlsdir);
if(ret<0) {
printf("[moldyn] visual init failure\n");
return ret;
/* calculate initial forces */
potential_force_calc(moldyn);
#ifdef DEBUG
-return 0;
+//return 0;
#endif
/* some stupid checks before we actually start calculating bullshit */
}
if(v) {
if(!(moldyn->total_steps%v)) {
- visual_atoms(&(moldyn->vis),moldyn->time,
- moldyn->atom,moldyn->count);
+ visual_atoms(moldyn);
}
}
/* display progress */
- if(!(moldyn->total_steps%10)) {
+ //if(!(moldyn->total_steps%10)) {
/* get current time */
gettimeofday(&t2,NULL);
/* copy over time */
t1=t2;
- }
+ //}
/* increase absolute time */
moldyn->time+=moldyn->tau;
}
#ifdef DEBUG
- printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
+ //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
+ if(moldyn->time>DSTART&&moldyn->time<DEND) {
+ printf("force:\n");
+ printf(" x: %0.40f\n",moldyn->atom[5832].f.x);
+ printf(" y: %0.40f\n",moldyn->atom[5832].f.y);
+ printf(" z: %0.40f\n",moldyn->atom[5832].f.z);
+ }
#endif
/* calculate global virial */
* restore function
*/
+int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
+
+ int fd;
+ int cnt,size;
+
+ fd=open(file,O_RDONLY);
+ if(fd<0) {
+ perror("[moldyn] load save file open");
+ return fd;
+ }
+
+ size=sizeof(t_moldyn);
+ cnt=read(fd,moldyn,size);
+ if(cnt!=size) {
+ perror("[moldyn] load save file read (moldyn)");
+ return cnt;
+ }
+
+ size=moldyn->count*sizeof(t_atom);
+
+ moldyn->atom=(t_atom *)malloc(size);
+ if(moldyn->atom==NULL) {
+ perror("[moldyn] load save file malloc (atoms)");
+ return -1;
+ }
+
+ cnt=read(fd,moldyn->atom,size);
+ if(cnt!=size) {
+ perror("[moldyn] load save file read (atoms)");
+ return cnt;
+ }
+
+ // hooks
+
+ return 0;
+}
+
int moldyn_load(t_moldyn *moldyn) {
// later ...
}
}
+int pair_correlation_init(t_moldyn *moldyn,double dr) {
+
+
+ return 0;
+}
+
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+ int slots;
+ double *stat;
+ int i,j;
+ t_linkcell *lc;
+ t_list neighbour[27];
+ t_atom *itom,*jtom;
+ t_list *this;
+ unsigned char bc;
+ t_3dvec dist;
+ double d,norm;
+ int o,s;
+ unsigned char ibrand;
+
+ lc=&(moldyn->lc);
+
+ slots=(int)(moldyn->cutoff/dr);
+ o=2*slots;
+
+ if(ptr!=NULL) {
+ stat=(double *)ptr;
+ }
+ else {
+ stat=(double *)malloc(3*slots*sizeof(double));
+ if(stat==NULL) {
+ perror("[moldyn] pair correlation malloc");
+ return -1;
+ }
+ }
+
+ memset(stat,0,3*slots*sizeof(double));
+
+ link_cell_init(moldyn,VERBOSE);
+
+ itom=moldyn->atom;
+
+ for(i=0;i<moldyn->count;i++) {
+ /* neighbour indexing */
+ link_cell_neighbour_index(moldyn,
+ (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (itom[i].r.y+moldyn->dim.y/2)/lc->x,
+ (itom[i].r.z+moldyn->dim.z/2)/lc->x,
+ neighbour);
+
+ /* brand of atom i */
+ ibrand=itom[i].brand;
+
+ for(j=0;j<27;j++) {
+ /* prepare the neighbour cell list */
+ this=&(neighbour[j]);
+ list_reset_f(this);
+
+ /* check for atoms */
+ if(this->start==NULL)
+ continue;
+
+ /* boundary check */
+ bc=(j<lc->dnlc)?0:1;
+
+ do {
+ jtom=this->current->data;
+
+
+ if(jtom==&(itom[i]))
+ continue;
+
+ /* only count pairs once */
+ if(itom[i].tag>jtom->tag)
+ continue;
+
+ /*
+ * pair correlation calc
+ */
+
+ /* distance */
+ v3_sub(&dist,&(jtom->r),&(itom[i].r));
+ if(bc) check_per_bound(moldyn,&dist);
+ d=v3_absolute_square(&dist);
+
+ /* ignore if greater cutoff */
+ if(d>moldyn->cutoff_square)
+ continue;
+
+ /* fill the slots */
+ d=sqrt(d);
+ s=(int)(d/dr);
+
+ if(ibrand!=jtom->brand) {
+ /* mixed */
+ stat[s]+=1;
+ }
+ else {
+ /* type a - type a bonds */
+ if(ibrand==0)
+ stat[s+slots]+=1;
+ else
+ /* type b - type b bonds */
+ stat[s+o]+=1;
+ }
+
+ } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+ }
+ }
+
+ /* normalization */
+ for(i=1;i<slots;i++) {
+ /*
+ * normalization: 4 pi r r dr
+ * here: not double counting pairs -> 2 pi r r dr
+ */
+ norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr;
+ stat[i]/=norm;
+ stat[slots+i]/=norm;
+ stat[o+i]/=norm;
+ }
+
+ if(ptr==NULL) {
+ /* todo: store/print pair correlation function */
+ free(stat);
+ }
+
+ free(moldyn->atom);
+
+ link_cell_shutdown(moldyn);
+
+ return 0;
+}
+
int analyze_bonds(t_moldyn *moldyn) {
* visualization code
*/
-int visual_init(t_visual *v,char *filebase) {
-
- char file[128+8];
+int visual_init(t_moldyn *moldyn,char *filebase) {
- strncpy(v->fb,filebase,128);
- memset(file,0,128+8);
+ strncpy(moldyn->vis.fb,filebase,128);
return 0;
}
-int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
+int visual_atoms(t_moldyn *moldyn) {
- int i,fd;
+ int i,j,fd;
char file[128+64];
t_3dvec dim;
double help;
+ t_visual *v;
+ t_atom *atom;
+ t_atom *btom;
+ t_linkcell *lc;
+ t_list neighbour[27];
+ u8 bc;
+ t_3dvec dist;
+ double d2;
+ u8 brand;
+ v=&(moldyn->vis);
dim.x=v->dim.x;
dim.y=v->dim.y;
dim.z=v->dim.z;
+ atom=moldyn->atom;
+ lc=&(moldyn->lc);
help=(dim.x+dim.y);
- sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,time);
+ sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
if(fd<0) {
perror("open visual save file fd");
}
/* write the actual data file */
+
+ // povray header
dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
- n,time,help/40.0,help/40.0,-0.8*help);
- for(i=0;i<n;i++)
+ moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
+
+ // atomic configuration
+ for(i=0;i<moldyn->count;i++) {
+ // atom type, positions, color and kinetic energy
dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
atom[i].r.x,
atom[i].r.y,
atom[i].r.z,
pse_col[atom[i].element],
atom[i].ekin);
+
+ /*
+ * bond detection should usually be done by potential
+ * functions. brrrrr! EVIL!
+ *
+ * todo: potentials need to export a 'find_bonds' function!
+ */
+
+ // bonds between atoms
+ if(!(atom[i].attr&ATOM_ATTR_VB))
+ continue;
+ link_cell_neighbour_index(moldyn,
+ (atom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (atom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (atom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour);
+ for(j=0;j<27;j++) {
+ list_reset_f(&neighbour[j]);
+ if(neighbour[j].start==NULL)
+ continue;
+ bc=j<lc->dnlc?0:1;
+ do {
+ btom=neighbour[j].current->data;
+ if(btom==&atom[i]) // skip identical atoms
+ continue;
+ //if(btom<&atom[i]) // skip half of them
+ // continue;
+ v3_sub(&dist,&(atom[i].r),&(btom->r));
+ if(bc) check_per_bound(moldyn,&dist);
+ d2=v3_absolute_square(&dist);
+ brand=atom[i].brand;
+ if(brand==btom->brand) {
+ if(d2>moldyn->bondlen[brand])
+ continue;
+ }
+ else {
+ if(d2>moldyn->bondlen[2])
+ continue;
+ }
+ dprintf(fd,"# [B] %f %f %f %f %f %f\n",
+ atom[i].r.x,atom[i].r.y,atom[i].r.z,
+ btom->r.x,btom->r.y,btom->r.z);
+ } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
+ }
+ }
+
+ // boundaries
if(dim.x) {
dprintf(fd,"# [D] %f %f %f %f %f %f\n",
-dim.x/2,-dim.y/2,-dim.z/2,