projects
/
physik
/
posic.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
create_lattice fixes, still virial probs ...
[physik/posic.git]
/
moldyn.c
diff --git
a/moldyn.c
b/moldyn.c
index
0d5027c
..
f6c3c81
100644
(file)
--- a/
moldyn.c
+++ b/
moldyn.c
@@
-339,22
+339,24
@@
int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
}
moldyn->atom=ptr;
atom=&(moldyn->atom[count]);
}
moldyn->atom=ptr;
atom=&(moldyn->atom[count]);
-
- v3_zero(&origin);
+
+ /* no atoms on the boundaries (only reason: it looks better!) */
+ origin.x=0.5*lc;
+ origin.y=0.5*lc;
+ origin.z=0.5*lc;
switch(type) {
case CUBIC:
set_nn_dist(moldyn,lc);
switch(type) {
case CUBIC:
set_nn_dist(moldyn,lc);
- origin.x=0.5*lc;
- origin.y=0.5*lc;
- origin.z=0.5*lc;
ret=cubic_init(a,b,c,lc,atom,&origin);
break;
case FCC:
ret=cubic_init(a,b,c,lc,atom,&origin);
break;
case FCC:
+ v3_scale(&origin,&origin,0.5);
set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
- ret=fcc_init(a,b,c,lc,atom,
NULL
);
+ ret=fcc_init(a,b,c,lc,atom,
&origin
);
break;
case DIAMOND:
break;
case DIAMOND:
+ v3_scale(&origin,&origin,0.25);
set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
ret=diamond_init(a,b,c,lc,atom,&origin);
break;
set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
ret=diamond_init(a,b,c,lc,atom,&origin);
break;
@@
-429,65
+431,55
@@
int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
int count;
int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
int count;
- int i,j;
+ int i,j
,k,l
;
t_3dvec o,r,n;
t_3dvec basis[3];
t_3dvec o,r,n;
t_3dvec basis[3];
- double help[3];
- double x,y,z;
-
- x=a*lc;
- y=b*lc;
- z=c*lc;
- if(origin) v3_copy(&o,origin);
- else v3_zero(&o);
+ count=0;
+ if(origin)
+ v3_copy(&o,origin);
+ else
+ v3_zero(&o);
/* construct the basis */
/* construct the basis */
- for(i=0;i<3;i++) {
- for(j=0;j<3;j++) {
-
if(i!=j) help[j]
=0.5*lc;
-
else help[j]=.0
;
- }
-
v3_set(&basis[i],help)
;
- }
+ memset(basis,0,3*sizeof(t_3dvec));
+ basis[0].x=0.5*lc;
+
basis[0].y
=0.5*lc;
+
basis[1].x=0.5*lc
;
+ basis[1].z=0.5*lc;
+
basis[2].y=0.5*lc
;
+ basis[2].z=0.5*lc;
- v3_zero(&r);
- count=0;
-
/* fill up the room */
r.x=o.x;
/* fill up the room */
r.x=o.x;
-
while(r.x<x
) {
+
for(i=0;i<a;i++
) {
r.y=o.y;
r.y=o.y;
-
while(r.y<y
) {
+
for(j=0;j<b;j++
) {
r.z=o.z;
r.z=o.z;
- while(r.z<z) {
+ for(k=0;k<c;k++) {
+ /* first atom */
v3_copy(&(atom[count].r),&r);
v3_copy(&(atom[count].r),&r);
- atom[count].element=1;
count+=1;
count+=1;
- for(i=0;i<3;i++) {
- v3_add(&n,&r,&basis[i]);
- if((n.x<x+o.x)&&
- (n.y<y+o.y)&&
- (n.z<z+o.z)) {
- v3_copy(&(atom[count].r),&n);
- count+=1;
- }
+ r.z+=lc;
+ /* the three face centered atoms */
+ for(l=0;l<3;l++) {
+ v3_add(&n,&r,&basis[l]);
+ v3_copy(&(atom[count].r),&n);
+ count+=1;
}
}
- r.z+=lc;
}
r.y+=lc;
}
r.x+=lc;
}
}
r.y+=lc;
}
r.x+=lc;
}
-
+
/* coordinate transformation */
/* coordinate transformation */
- help[0]=x/2.0;
- help[1]=y/2.0;
- help[2]=z/2.0;
- v3_set(&n,help);
- for(i=0;i<count;i++)
- v3_sub(&(atom[i].r),&(atom[i].r),&n);
-
+ for(i=0;i<count;i++) {
+ atom[i].r.x-=(a*lc)/2.0;
+ atom[i].r.y-=(b*lc)/2.0;
+ atom[i].r.z-=(c*lc)/2.0;
+ }
+
return count;
}
return count;
}
@@
-1575,6
+1567,7
@@
int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
v3_scale(&force,&force,-1.0); /* f = - grad E */
v3_add(&(ai->f),&(ai->f),&force);
virial_calc(ai,&force,&distance);
v3_scale(&force,&force,-1.0); /* f = - grad E */
v3_add(&(ai->f),&(ai->f),&force);
virial_calc(ai,&force,&distance);
+if(force.x*distance.x<=0) printf("virial xx: %.15f -> %f %f %f\n",force.x*distance.x,distance.x,distance.y,distance.z);
virial_calc(aj,&force,&distance); /* f and d signe switched */
}
virial_calc(aj,&force,&distance); /* f and d signe switched */
}