pthread_mutex_t emutex;
#endif
+/* fully constrained relaxation technique - global pointers */
+u8 crtt;
+u8 *constraints;
+double *trafo_angle;
+
/*
* the moldyn functions
*/
pthread_mutex_init(&emutex,NULL);
#endif
-#ifdef CONSTRAINT_110_5832
- printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n");
- printf("\n\n\n!! -- constraints enabled -- !!\n\n\n");
-#endif
-#ifdef CONSTRAINT_11X_5832
- printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n");
- printf("\n\n\n!! -- constraints enabled -- !!\n\n\n");
-#endif
return 0;
}
return 0;
}
+/* basis trafo */
+
+#define FORWARD 0
+#define BACKWARD 1
+
+int basis_trafo(t_3dvec *r,u8 dir,double z,double x) {
+
+ t_3dvec tmp;
+
+ if(dir==FORWARD) {
+ if(z!=0.0) {
+ v3_copy(&tmp,r);
+ r->x=cos(z)*tmp.x-sin(z)*tmp.y;
+ r->y=sin(z)*tmp.x+cos(z)*tmp.y;
+ }
+ if(x!=0.0) {
+ v3_copy(&tmp,r);
+ r->y=cos(x)*tmp.y-sin(x)*tmp.z;
+ r->z=sin(x)*tmp.y+cos(x)*tmp.z;
+ }
+ }
+ else {
+ if(x!=0.0) {
+ v3_copy(&tmp,r);
+ r->y=cos(-x)*tmp.y-sin(-x)*tmp.z;
+ r->z=sin(-x)*tmp.y+cos(-x)*tmp.z;
+ }
+ if(z!=0.0) {
+ v3_copy(&tmp,r);
+ r->x=cos(-z)*tmp.x-sin(-z)*tmp.y;
+ r->y=sin(-z)*tmp.x+cos(-z)*tmp.y;
+ }
+ }
+
+ return 0;
+}
+
/* velocity verlet */
int velocity_verlet(t_moldyn *moldyn) {
double tau,tau_square,h;
t_3dvec delta;
t_atom *atom;
-#ifdef CONSTRAINT_11X_5832
- double xt,yt,zt;
- double xtt,ytt,ztt;
-#endif
atom=moldyn->atom;
count=moldyn->count;
tau=moldyn->tau;
tau_square=moldyn->tau_square;
-#ifdef CONSTRAINT_110_5832
- if(count==5833) {
- atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
- atom[5832].f.y=-atom[5832].f.x;
- }
-#endif
-#ifdef CONSTRAINT_11X_5832
- if(count==5833) {
- // second trafo
- xt=atom[5832].f.x;
- yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129);
- zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129);
- // first trafo
- xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
- ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
- ztt=zt;
- // apply constraints
- ytt=0.0;
- // first trafo backwards
- xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
- yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
- zt=ztt;
- // second trafo backwards
- atom[5832].f.x=xt;
- atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
- atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
- }
-#endif
for(i=0;i<count;i++) {
+
/* check whether fixed atom */
if(atom[i].attr&ATOM_ATTR_FP)
continue;
+
/* new positions */
h=0.5/atom[i].mass;
- v3_scale(&delta,&(atom[i].v),tau);
-#ifdef CONSTRAINT_110_5832
- if(i==5832) {
- delta.y=-delta.x;
+
+ /* constraint relaxation */
+ if(crtt) {
+ // forces
+ basis_trafo(&(atom[i].f),FORWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ if(constraints[3*i])
+ atom[i].f.x=0;
+ if(constraints[3*i+1])
+ atom[i].f.y=0;
+ if(constraints[3*i+2])
+ 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]);
}
-#endif
-#ifdef JHDVHJDV
-#ifdef CONSTRAINT_11X_5832
- if(i==5833) {
- // second trafo
- xt=delta.x;
- yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
- zt=delta.y*sin(-0.16935129)+delta.z*cos(-0.16935129);
- // first trafo
- xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
- ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
- ztt=zt;
- // apply constraints
- ytt=0.0;
- // first trafo backwards
- xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
- yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
- zt=ztt;
- // second trafo backwards
- delta.x=xt;
- delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
- delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
- }
-#endif
-#endif
+
#ifndef QUENCH
+ v3_scale(&delta,&(atom[i].v),tau);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
#endif
v3_scale(&delta,&(atom[i].f),h*tau_square);
-#ifdef CONSTRAINT_110_5832
- if(i==5832) {
- delta.y=-delta.x;
- }
-#endif
-#ifdef GHDJHBSJBSD
-#ifdef CONSTRAINT_11X_5832
- if(i==5833) {
- // second trafo
- xt=delta.x;
- yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
- zt=delta.y*sin(-0.16935129)+delta.z*cos(-0.16935129);
- // first trafo
- xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
- ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
- ztt=zt;
- // apply constraints
- ytt=0.0;
- // first trafo backwards
- xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
- yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
- zt=ztt;
- // second trafo backwards
- delta.x=xt;
- delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
- delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
- }
-#endif
-#endif
v3_add(&(atom[i].r),&(atom[i].r),&delta);
//check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
check_per_bound(moldyn,&(atom[i].r));
albe_potential_force_calc(moldyn);
#endif
-#ifdef CONSTRAINT_110_5832
- if(count==5833) {
- atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
- atom[5832].f.y=-atom[5832].f.x;
- }
-#endif
-#ifdef CONSTRAINT_11X_5832
- if(count==5833) {
- // second trafo
- xt=atom[5832].f.x;
- yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129);
- zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129);
- // first trafo
- xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
- ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
- ztt=zt;
- // apply constraints
- ytt=0.0;
- // first trafo backwards
- xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
- yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
- zt=ztt;
- // second trafo backwards
- atom[5832].f.x=xt;
- atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
- atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
- }
-#endif
for(i=0;i<count;i++) {
/* check whether fixed atom */
if(atom[i].attr&ATOM_ATTR_FP)
continue;
+
+ /* constraint relaxation */
+ if(crtt) {
+ // forces
+ basis_trafo(&(atom[i].f),FORWARD,
+ trafo_angle[2*i],trafo_angle[2*i+1]);
+ if(constraints[3*i])
+ atom[i].f.x=0;
+ if(constraints[3*i+1])
+ atom[i].f.y=0;
+ if(constraints[3*i+2])
+ 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)] */
v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
v3_add(&(atom[i].v),&(atom[i].v),&delta);