}
int set_temperature(t_moldyn *moldyn,double t) {
-
+
moldyn->t=t;
return 0;
int i;
double e,scale;
t_atom *atom;
+ int count;
atom=moldyn->atom;
*/
e=0.0;
- for(i=0;i<moldyn->count;i++)
- e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
- scale=(1.5*moldyn->count*K_BOLTZMANN*moldyn->t)/e;
+ count=0;
+ for(i=0;i<moldyn->count;i++) {
+ if(atom[i].attr&ATOM_ATTR_HB) {
+ e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+ count+=1;
+ }
+ }
+
+ /* temporary hack for e,t = 0 */
+ if(e==0.0) {
+ if(moldyn->t!=0.0)
+ thermal_init(moldyn);
+ else
+ return 0;
+ }
+
+ /* direct scaling */
+ scale=(1.5*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),scale);
+ if(atom[i].attr&ATOM_ATTR_HB)
+ v3_scale(&(atom[i].v),&(atom[i].v),scale);
return 0;
}
count2=27;
a=nx*ny;
-
cell[0]=lc->subcell[i+j*nx+k*a];
for(ci=-1;ci<=1;ci++) {
bx=0;
/* zero absolute time */
moldyn->time=0.0;
-
for(sched=0;sched<moldyn->schedule.content_count;sched++) {
/* setting amount of runs and finite time step size */
int potential_force_calc(t_moldyn *moldyn) {
int i,j,k,count;
- t_atom *atom,*btom,*ktom;
+ t_atom *itom,*jtom,*ktom;
t_linkcell *lc;
- t_list neighbour[27];
- t_list *this,*thisk,*neighbourk;
- u8 bc,bck;
+ t_list neighbour_i[27],neighbour_j[27];
+ t_list *this,*that;
+ u8 bc_ij,bc_ijk;
int countn,dnlc;
count=moldyn->count;
- atom=moldyn->atom;
+ itom=moldyn->atom;
lc=&(moldyn->lc);
/* reset energy */
for(i=0;i<count;i++) {
/* reset force */
- v3_zero(&(atom[i].f));
+ v3_zero(&(itom[i].f));
/* single particle potential/force */
- if(atom[i].attr&ATOM_ATTR_1BP)
- moldyn->func1b(moldyn,&(atom[i]));
+ if(itom[i].attr&ATOM_ATTR_1BP)
+ moldyn->func1b(moldyn,&(itom[i]));
/* 2 body pair potential/force */
- if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+ if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
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);
+ (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (itom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (itom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour_i);
countn=lc->countn;
dnlc=lc->dnlc;
for(j=0;j<countn;j++) {
- this=&(neighbour[j]);
+ this=&(neighbour_i[j]);
list_reset(this);
if(this->start==NULL)
continue;
- bc=(j<dnlc)?0:1;
+ bc_ij=(j<dnlc)?0:1;
do {
- btom=this->current->data;
+ jtom=this->current->data;
- if(btom==&(atom[i]))
+ if(jtom==&(itom[i]))
continue;
- if((btom->attr&ATOM_ATTR_2BP)&
- (atom[i].attr&ATOM_ATTR_2BP))
+ if((jtom->attr&ATOM_ATTR_2BP)&
+ (itom[i].attr&ATOM_ATTR_2BP))
moldyn->func2b(moldyn,
- &(atom[i]),
- btom,
- bc);
+ &(itom[i]),
+ jtom,
+ bc_ij);
/* 3 body potential/force */
- if(!(atom[i].attr&ATOM_ATTR_3BP)||
- !(btom->attr&ATOM_ATTR_3BP))
+ if(!(itom[i].attr&ATOM_ATTR_3BP)||
+ !(jtom->attr&ATOM_ATTR_3BP))
continue;
-printf("DEBUG: problem exists here ...\n");
link_cell_neighbour_index(moldyn,
- (btom->r.x+moldyn->dim.x/2)/lc->x,
- (btom->r.y+moldyn->dim.y/2)/lc->y,
- (btom->r.z+moldyn->dim.z/2)/lc->z,
- neighbourk);
-printf("DEBUG: as you won't see that!\n");
+ (jtom->r.x+moldyn->dim.x/2)/lc->x,
+ (jtom->r.y+moldyn->dim.y/2)/lc->y,
+ (jtom->r.z+moldyn->dim.z/2)/lc->z,
+ neighbour_j);
+ /* neighbours of j */
for(k=0;k<lc->countn;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
+ that=&(neighbour_j[k]);
+ list_reset(that);
+
+ if(that->start==NULL)
+ continue;
+
+ bc_ijk=(k<lc->dnlc)?0:1;
+
+ do {
+
+ ktom=that->current->data;
+
+ if(!(ktom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ if(ktom==jtom)
+ continue;
+
+ if(ktom==&(itom[i]))
+ continue;
+
+ moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+
+ } while(list_next(that)!=\
+ L_NO_NEXT_ELEMENT);
+
+ }
+
+ /* neighbours of i */
+ for(k=0;k<countn;k++) {
+
+ that=&(neighbour_i[k]);
+ list_reset(that);
- if(thisk->start==NULL)
+ if(that->start==NULL)
continue;
- bck=(k<lc->dnlc)?0:1;
+ bc_ijk=(k<dnlc)?0:1;
do {
- ktom=thisk->current->data;
+ ktom=that->current->data;
if(!(ktom->attr&ATOM_ATTR_3BP))
continue;
- if(ktom==btom)
+ if(ktom==jtom)
continue;
- if(ktom==&(atom[i]))
+ if(ktom==&(itom[i]))
continue;
- moldyn->func3b(moldyn,&(atom[i]),btom,ktom,bck);
+ moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
- } while(list_next(thisk)!=\
+ } while(list_next(that)!=\
L_NO_NEXT_ELEMENT);
}
projects / physik / posic.git / blobdiff