lots of mostly small changes ...

diff --git a/Makefile b/Makefile
index e171354..5f25159 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -8,5 +8,6 @@ all: posic
 posic: $(OBJS) moldyn.o
        $(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c
 
+.PHONY:clean
 clean:
        rm -f *.o posic */*.o

diff --git a/moldyn.c b/moldyn.c
index 0fd1aca..db575e9 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -547,7 +547,7 @@ double potential_lennard_jones(t_moldyn *moldyn) {
        params=moldyn->pot_params;
        atom=moldyn->atom;
        count=moldyn->count;
-       eps=params->epsilon;
+       eps=params->epsilon4;
        sig6=params->sigma6;
        sig12=params->sigma12;
 
@@ -559,7 +559,7 @@ double potential_lennard_jones(t_moldyn *moldyn) {
                        help=d*d;                               /* 1/r^4 */
                        help*=d;                                /* 1/r^6 */
                        d=help*help;                            /* 1/r^12 */
-                       u+=eps*(sig12*d-sig6*help);
+                       u+=eps*(sig6*help-sig12*d);
                }
        }
        
@@ -579,7 +579,7 @@ int force_lennard_jones(t_moldyn *moldyn) {
        atom=moldyn->atom;      
        count=moldyn->count;
        params=moldyn->pot_params;
-       eps=params->epsilon;
+       eps=params->epsilon4;
        sig6=params->sigma6;
        sig12=params->sigma12;
 
@@ -598,11 +598,13 @@ int force_lennard_jones(t_moldyn *moldyn) {
                                h1*=h2;                         /* 1/r^14 */
                                h1*=sig12;
                                h2*=sig6;
-                               d=12.0*h1-6.0*h2;
+                               /* actually there would be a '-',       *
+                                * but f=-d/dr potential                */
+                               d=-12.0*h1+6.0*h2;
                                d*=eps;
                                v3_scale(&force,&distance,d);
-                               v3_add(&(atom[j].f),&(atom[j].f),&force);
-                               v3_sub(&(atom[i].f),&(atom[i].f),&force);
+                               v3_add(&(atom[i].f),&(atom[i].f),&force);
+                               v3_sub(&(atom[j].f),&(atom[j].f),&force);
                        }
                }
        }

diff --git a/moldyn.h b/moldyn.h
index 42764c5..266436b 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -22,6 +22,8 @@ typedef struct s_atom {
        //t_list vicinity       /* verlet neighbour list */
 } t_atom;
 
+#include "visual/visual.h"
+
 typedef struct s_moldyn {
        /* atoms, amount, dimensions */
        int count;
@@ -68,7 +70,7 @@ typedef struct s_ho_params {
 typedef struct s_lj_params {
        double sigma6;
        double sigma12;
-       double epsilon;
+       double epsilon4;
 } t_lj_params;
 
 /*

diff --git a/posic.c b/posic.c
index f026abe..992f9d1 100644 (file)
--- a/posic.c
+++ b/posic.c
@@ -53,14 +53,14 @@ int main(int argc,char **argv) {
        vis.dim.y=b*LC_SI;
        vis.dim.z=c*LC_SI;
 
-       /* init lattice
+       /* init lattice */
        printf("placing silicon atoms ... ");
        count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
-       printf("(%d) ok!\n",count); */
-       /* testing purpose */
+       printf("(%d) ok!\n",count);
+       /* testing purpose
        count=2;
        si=malloc(2*sizeof(t_atom));
-       si[0].r.x=0.35*sqrt(3.0)*LC_SI/2.0;
+       si[0].r.x=0.23*sqrt(3.0)*LC_SI/2.0;
        si[0].r.y=0;
        si[0].r.z=0;
        si[0].element=SI;
@@ -70,7 +70,7 @@ int main(int argc,char **argv) {
        si[1].r.z=0;
        si[1].element=SI;
        si[1].mass=M_SI;
-       /* */
+       */
 
        /* moldyn init (now si is a valid address) */
        md.count=count;
@@ -90,10 +90,8 @@ int main(int argc,char **argv) {
        md.visual=&vis;
 
        printf("setting thermal fluctuations (T=%f K)\n",md.t);
-       //thermal_init(&md,&random,count);
-       for(a=0;a<count;a++) v3_zero(&(si[0].v));
-       //v3_zero(&(si[0].v));
-       //v3_zero(&(si[1].v));
+       thermal_init(&md,&random,count);
+       //for(a=0;a<count;a++) v3_zero(&(si[0].v));
 
        /* check kinetic energy */
 
@@ -110,10 +108,10 @@ int main(int argc,char **argv) {
        help=lj.sigma6*lj.sigma6;
        lj.sigma6*=help;
        lj.sigma12=lj.sigma6*lj.sigma6;
-       lj.epsilon=LJ_EPSILON_SI;
+       lj.epsilon4=4.0*LJ_EPSILON_SI;
 
        ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
-       ho.spring_constant=LJ_EPSILON_SI;
+       ho.spring_constant=4.0*LJ_EPSILON_SI;
 
        u=get_e_pot(&md);
 

diff --git a/run b/run
new file mode 100755 (executable)
index 0000000..4c753fe
--- /dev/null
+++ b/run
@@ -0,0 +1,12 @@
+./clean
+./posic $@
+if [ "$?" == "0" ]; then
+       ./perms
+       argv="`echo $@ | grep -- '-V'`"
+       if [ "$argv" ] ; then
+               base=`echo $argv | awk -F'-V' '{ print $2 }' | \
+                       awk '{ print $2 }'`
+               rasmol -nodisplay < ${base}.scr
+               ./ppm2avi
+       fi
+fi

diff --git a/visual/visual.c b/visual/visual.c
index 38ff0b0..29d9469 100644 (file)
--- a/visual/visual.c
+++ b/visual/visual.c
@@ -17,6 +17,29 @@
 #include "../moldyn.h"
 #include "../math/math.h"
 
+/* the pse, only needed here */
+static char *pse[]={
+       "*",
+       "H",
+       "He",
+       "Li",
+       "Be",
+       "B",
+       "C",
+       "N",
+       "O",
+       "F",
+       "Ne",
+       "Na",
+       "Mg",
+       "Al",
+       "Si",
+       "P",
+       "S",
+       "Cl",
+       "Ar",
+};
+
 int visual_init(t_visual *v,char *filebase) {
 
        char file[128+8];
@@ -30,6 +53,7 @@ int visual_init(t_visual *v,char *filebase) {
                perror("open visual fd");
                return -1;
        }
+       dprintf(v->fd,"set write on\n");
 
        memset(&(v->dim),0,sizeof(t_3dvec));
 
@@ -38,7 +62,10 @@ int visual_init(t_visual *v,char *filebase) {
 
 int visual_tini(t_visual *v) {
 
-       if(v->fd) close(v->fd);
+       if(v->fd) {
+               dprintf(v->fd,"quit\n");
+               close(v->fd);
+       }
 
        return 0;
 }
@@ -53,28 +80,6 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        dim.y=10e9*v->dim.y;
        dim.z=10e9*v->dim.z;
 
-       char pse[19][4]={
-               "*",
-               "H",
-               "He",
-               "Li",
-               "Be",
-               "B",
-               "C",
-               "N",
-               "O",
-               "F",
-               "Ne",
-               "Na",
-               "Mg",
-               "Al",
-               "Si",
-               "P",
-               "S",
-               "Cl",
-               "Ar",
-       };
-
        sprintf(file,"%s-%.15f.xyz",v->fb,time);
        fd=open(file,O_WRONLY|O_CREAT|O_TRUNC);
        if(fd<0) {
@@ -84,7 +89,7 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
 
        /* script file update */
        dprintf(v->fd,"load xyz %s\n",file);
-       dprintf(v->fd,"spacefill 300\n");
+       dprintf(v->fd,"spacefill 270\n");
        dprintf(v->fd,"rotate x 100\n");
        dprintf(v->fd,"rotate y 10\n");
        dprintf(v->fd,"set ambient 20\n");

diff --git a/visual/visual.h b/visual/visual.h
index e2681d9..8726741 100644 (file)
--- a/visual/visual.h
+++ b/visual/visual.h
@@ -11,13 +11,14 @@
 #include "../moldyn.h"
 #include "../math/math.h"
 
+/* datatypes */
+
 typedef struct s_visual {
        int fd;                 /* rasmol script file descriptor */
        char fb[128];           /* basename of the save files */
        t_3dvec dim;            /* dimensions of the simulation cell */
 } t_visual;
 
-
 /* function prototypes */
 
 int visual_init(t_visual *v,char *filebase);
@@ -25,3 +26,5 @@ int visual_tini(t_visual *v);
 int visual_atoms(t_visual *v,double time,t_atom *atom,int n);
 
 #endif
+
+