added Verlet scheme and NxN non-bonded functionality

[alexxy/gromacs.git] / src / mdlib / shellfc.c

diff --git a/src/mdlib/shellfc.c b/src/mdlib/shellfc.c
index 5d9e7febe06264253742ee3c5860ee23d5a18dc0..86ee3c05aa58c64890beeda8b72fda638d5afe71 100644 (file)
--- a/src/mdlib/shellfc.c
+++ b/src/mdlib/shellfc.c
@@ -117,7 +117,12 @@ static void predict_shells(FILE *fplog,rvec x[],rvec v[],real dt,
   int  i,m,s1,n1,n2,n3;
   real dt_1,dt_2,dt_3,fudge,tm,m1,m2,m3;
   rvec *ptr;
+  gmx_mtop_atomlookup_t alook=NULL;
   t_atom *atom;
+
+  if (mass == NULL) {
+    alook = gmx_mtop_atomlookup_init(mtop);
+  }
   
   /* We introduce a fudge factor for performance reasons: with this choice
    * the initial force on the shells is about a factor of two lower than 
@@ -171,11 +176,11 @@ static void predict_shells(FILE *fplog,rvec x[],rvec v[],real dt,
        m3 = mass[n3];
       } else {
        /* Not the correct masses with FE, but it is just a prediction... */
-       gmx_mtop_atomnr_to_atom(mtop,n1,&atom);
+       gmx_mtop_atomnr_to_atom(alook,n1,&atom);
        m1 = atom->m;
-       gmx_mtop_atomnr_to_atom(mtop,n2,&atom);
+       gmx_mtop_atomnr_to_atom(alook,n2,&atom);
        m2 = atom->m;
-       gmx_mtop_atomnr_to_atom(mtop,n3,&atom);
+       gmx_mtop_atomnr_to_atom(alook,n3,&atom);
        m3 = atom->m;
       }
       tm = dt_1/(m1+m2+m3);
@@ -186,6 +191,10 @@ static void predict_shells(FILE *fplog,rvec x[],rvec v[],real dt,
       gmx_fatal(FARGS,"Shell %d has %d nuclei!",i,s[i].nnucl);
     }
   }
+
+  if (mass == NULL) {
+    gmx_mtop_atomlookup_destroy(alook);
+  }
 }
 
 gmx_shellfc_t init_shell_flexcon(FILE *fplog,
@@ -698,7 +707,8 @@ static void init_adir(FILE *log,gmx_shellfc_t shfc,
                      t_commrec *cr,int dd_ac1,
                      gmx_large_int_t step,t_mdatoms *md,int start,int end,
                      rvec *x_old,rvec *x_init,rvec *x,
-                     rvec *f,rvec *acc_dir,matrix box,
+                     rvec *f,rvec *acc_dir,
+                     gmx_bool bMolPBC,matrix box,
                      real *lambda,real *dvdlambda,t_nrnb *nrnb)
 {
   rvec   *xnold,*xnew;
@@ -740,13 +750,14 @@ static void init_adir(FILE *log,gmx_shellfc_t shfc,
     }
   }
   constrain(log,FALSE,FALSE,constr,idef,ir,NULL,cr,step,0,md,
-           x,xnold-start,NULL,box,
-           lambda[efptBONDED],&(dvdlambda[efptBONDED]),NULL,NULL,nrnb,econqCoord,FALSE,0,0);
+           x,xnold-start,NULL,bMolPBC,box,
+           lambda[efptBONDED],&(dvdlambda[efptBONDED]),
+           NULL,NULL,nrnb,econqCoord,FALSE,0,0);
   constrain(log,FALSE,FALSE,constr,idef,ir,NULL,cr,step,0,md,
-           x,xnew-start,NULL,box,
-           lambda[efptBONDED],&(dvdlambda[efptBONDED]),NULL,NULL,nrnb,econqCoord,FALSE,0,0);
+           x,xnew-start,NULL,bMolPBC,box,
+           lambda[efptBONDED],&(dvdlambda[efptBONDED]),
+           NULL,NULL,nrnb,econqCoord,FALSE,0,0);
 
-  /* Set xnew to minus the acceleration */
   for (n=start; n<end; n++) {
     for(d=0; d<DIM; d++)
       xnew[n-start][d] =
@@ -757,8 +768,9 @@ static void init_adir(FILE *log,gmx_shellfc_t shfc,
 
   /* Project the acceleration on the old bond directions */
   constrain(log,FALSE,FALSE,constr,idef,ir,NULL,cr,step,0,md,
-           x_old,xnew-start,acc_dir,box,
-           lambda[efptBONDED],&(dvdlambda[efptBONDED]),NULL,NULL,nrnb,econqDeriv_FlexCon,FALSE,0,0);
+           x_old,xnew-start,acc_dir,bMolPBC,box,
+           lambda[efptBONDED],&(dvdlambda[efptBONDED]),
+           NULL,NULL,nrnb,econqDeriv_FlexCon,FALSE,0,0); 
 }
 
 int relax_shell_flexcon(FILE *fplog,t_commrec *cr,gmx_bool bVerbose,
@@ -897,7 +909,8 @@ int relax_shell_flexcon(FILE *fplog,t_commrec *cr,gmx_bool bVerbose,
     init_adir(fplog,shfc,
              constr,idef,inputrec,cr,dd_ac1,mdstep,md,start,end,
              shfc->x_old-start,state->x,state->x,force[Min],
-             shfc->acc_dir-start,state->box,state->lambda,&dum,nrnb);
+             shfc->acc_dir-start,
+             fr->bMolPBC,state->box,state->lambda,&dum,nrnb);
 
     for(i=start; i<end; i++)
       sf_dir += md->massT[i]*norm2(shfc->acc_dir[i-start]);
@@ -952,7 +965,7 @@ int relax_shell_flexcon(FILE *fplog,t_commrec *cr,gmx_bool bVerbose,
       init_adir(fplog,shfc,
                constr,idef,inputrec,cr,dd_ac1,mdstep,md,start,end,
                x_old-start,state->x,pos[Min],force[Min],acc_dir-start,
-               state->box,state->lambda,&dum,nrnb);
+               fr->bMolPBC,state->box,state->lambda,&dum,nrnb);
       
       directional_sd(fplog,pos[Min],pos[Try],acc_dir-start,start,end,
                     fr->fc_stepsize);
@@ -986,7 +999,7 @@ int relax_shell_flexcon(FILE *fplog,t_commrec *cr,gmx_bool bVerbose,
       init_adir(fplog,shfc,
                constr,idef,inputrec,cr,dd_ac1,mdstep,md,start,end,
                x_old-start,state->x,pos[Try],force[Try],acc_dir-start,
-               state->box,state->lambda,&dum,nrnb);
+               fr->bMolPBC,state->box,state->lambda,&dum,nrnb);
 
       for(i=start; i<end; i++)
        sf_dir += md->massT[i]*norm2(acc_dir[i-start]);