oups, forgot the basis trafo stuff

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index d715c52..33224a6 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -2319,6 +2319,7 @@ int velocity_verlet(t_moldyn *moldyn) {
 
                /* constraint relaxation */
                if(crtt) {
 
                /* constraint relaxation */
                if(crtt) {

+                       // forces

                        basis_trafo(&(atom[i].f),FORWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);
                        if(constraints[3*i])
                        basis_trafo(&(atom[i].f),FORWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);
                        if(constraints[3*i])


@@ -2329,6 +2330,17 @@ int velocity_verlet(t_moldyn *moldyn) {
                                atom[i].f.z=0;
                        basis_trafo(&(atom[i].f),BACKWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);
                                atom[i].f.z=0;
                        basis_trafo(&(atom[i].f),BACKWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);

+                       // velocities
+                       basis_trafo(&(atom[i].v),FORWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+                       if(constraints[3*i])
+                               atom[i].v.x=0;
+                       if(constraints[3*i+1])
+                               atom[i].v.y=0;
+                       if(constraints[3*i+2])
+                               atom[i].v.z=0;
+                       basis_trafo(&(atom[i].v),BACKWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);

                }
 
 #ifndef QUENCH
                }
 
 #ifndef QUENCH


@@ -2369,6 +2381,7 @@ int velocity_verlet(t_moldyn *moldyn) {
 
                /* constraint relaxation */
                if(crtt) {
 
                /* constraint relaxation */
                if(crtt) {

+                       // forces

                        basis_trafo(&(atom[i].f),FORWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);
                        if(constraints[3*i])
                        basis_trafo(&(atom[i].f),FORWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);
                        if(constraints[3*i])


@@ -2379,6 +2392,17 @@ int velocity_verlet(t_moldyn *moldyn) {
                                atom[i].f.z=0;
                        basis_trafo(&(atom[i].f),BACKWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);
                                atom[i].f.z=0;
                        basis_trafo(&(atom[i].f),BACKWARD,
                                    trafo_angle[2*i],trafo_angle[2*i+1]);

+                       // velocities
+                       basis_trafo(&(atom[i].v),FORWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);
+                       if(constraints[3*i])
+                               atom[i].v.x=0;
+                       if(constraints[3*i+1])
+                               atom[i].v.y=0;
+                       if(constraints[3*i+2])
+                               atom[i].v.z=0;
+                       basis_trafo(&(atom[i].v),BACKWARD,
+                                   trafo_angle[2*i],trafo_angle[2*i+1]);

                }
 
                /* again velocities [actually v(t+tau)] */
                }
 
                /* again velocities [actually v(t+tau)] */