--- trunk/ripple/dp.c	2004/05/24 16:42:34	1190
+++ trunk/ripple/dp.c	2004/05/24 16:52:32	1191
@@ -35,7 +35,9 @@ struct system {
   double z0;               // default z axis position
   double theta0;           // default theta angle
   double kz;               // force constant for z displacement
+  double kr;               // force constant for z displacement
   double ktheta;           // force constant for theta displacement
+  double t;                // the temperature of the system
   int    nCycles;          // How many cycles to do in total
   int    iCycle;           // How many cycles have we done?
   int    nMoves;           // How many MC moves in each cycle
@@ -50,6 +52,7 @@ struct system {
   int    nAccepts;         // number of MC moves that have been accepted
   int    nx;               // number of unit cells in x direction
   int    ny;               // number of unit cells in y direction
+  double XYNNDIST;         // maximum distance of nearest neighbors in XY plane
 };
 
 char *program_name; /* the name of the program */
@@ -129,12 +132,14 @@ int main(argc, argv) 
   state->strength = 7.0;
   state->z0 = 0.0;
   state->kz = kb; 
+  state->kr = kb;
   state->theta0 = M_PI / 2.0;
   state->ktheta = 0.0; 
   state->dtheta = 0.3;
   state->deltaz = 1.0;
   state->deltaphi = 0.5; 
-  state->beta = 1.0 / (kb * 300.0);
+  state->t = 300;
+  state->beta = 1.0 / (kb * state->t);
 
   // these three should really be given as command line arguments, but we'll
   // set defaults here just in case:
@@ -232,9 +237,25 @@ int main(argc, argv) 
 	  // -k kz   sets the value of kz in units of kb
 	  i++;
 	  state->kz = kb * atof( argv[i] );
+	  state->kr = state->kz;
 	  done = 1;
 	  break;
           
+	case 'q':
+	  // -q mu set the strength of the dipole
+	  i++;
+	  state->strength =  atof( argv[i] );
+	  done = 1;
+	  break;
+ 
+	case 't':
+	  // -t set the temperature of the system
+	  i++;
+	  state->t =  atof( argv[i] );
+	  state->beta = 1.0 / (kb * state->t);
+	  done = 1;
+	  break;
+
         case 'h':
           // -h hexSpace  set the inter-atomic spacing for a regular
           //              hexagonal lattice
@@ -391,6 +412,14 @@ int main(argc, argv) 
       aLat = dx / (double)(state->nx);
       bLat = dy / (double)(state->ny);
       
+      if (0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)) > bLat) {
+         state->XYNNDIST =  0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0));
+      } else {
+         state->XYNNDIST = bLat;
+      }
+      // slightly larger so we can use <  as a comparison 
+      state->XYNNDIST = state->XYNNDIST * 1.01;
+
       // Find HmatI:
       
       invertMat2(state->Hmat, state->HmatI);
@@ -509,6 +538,14 @@ int main(argc, argv) 
     state->nx = cells;
     state->ny = (int) ((double)cells * aLat / bLat );
 
+    if (0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0)) > bLat) {
+       state->XYNNDIST =  0.5*sqrt(pow(aLat,2.0) + pow(bLat,2.0));
+    } else {
+       state->XYNNDIST = bLat;
+    }
+    // slightly larger so we can use <  as a comparison 
+    state->XYNNDIST = state->XYNNDIST * 1.01;
+
     // each cell has 2 atoms (one at (0,0) and one at (a/2 , b/2))
     state->nAtoms = 2 * state->nx * state->ny;
     
@@ -576,12 +613,12 @@ int main(argc, argv) 
     exit(8);
   }
   
-  if (state->iCycle >= state->nCycles) {
+/*  if (state->iCycle >= state->nCycles) {
     printf("This configuration has already gone past the requested number\n"
            "of MC cycles!  Use the -n flag to request more MC cycles!\n");
     exit(8);
   }
-
+*/
   // Do the MC simulation (finally!)
   
   en = toterg(state);
@@ -701,15 +738,24 @@ double eneri(struct system *state, struct coords iTemp
 
 double eneri(struct system *state, struct coords iTemp, int i, int jb) {
   
-  double uxi, uyi, uzi, r, r2, r3, r5, uxj, uyj, uzj, rcut;
+  double uxi, uyi, uzi, rxy, rij, r, r2, r3, r5, uxj, uyj, uzj, rcut;
   double udotu, rdotui, rdotuj, vij, pre, vint;
+  double dx, dy, aLat, bLat;
   double eni, dz;
   double d[2];
   int j;
+
+  vint = 0.0;
   
   pre = 14.38362;
   
   rcut = 30.0;
+
+  dx = sqrt(pow(state->Hmat[0][0], 2) + pow(state->Hmat[1][0], 2));
+  dy = sqrt(pow(state->Hmat[0][1], 2) + pow(state->Hmat[1][1], 2));
+
+  aLat = dx / (double)(state->nx);
+  bLat = dy / (double)(state->ny);
   
   eni = 0.0;
   uxi = iTemp.mu * cos(iTemp.phi) * sin(iTemp.theta);
@@ -729,10 +775,16 @@ double eneri(struct system *state, struct coords iTemp
       // z is unwrapped!
 
       dz = state->r[j].pos[2] - iTemp.pos[2];
-          
+
+
+      rxy = sqrt(pow(d[0], 2) + pow(d[1], 2));
       r2 = pow(d[0], 2) + pow(d[1], 2) + pow(dz, 2);
       r = sqrt(r2);
 
+      if (rxy < state->XYNNDIST) {
+	 vint += 0.5 * state->kr * r2;
+      }
+
       if(r < rcut) {
 
         r3 = r2*r;
@@ -751,7 +803,6 @@ double eneri(struct system *state, struct coords iTemp
     } 
   }
   
-  vint = 0.5 * state->kz * pow((iTemp.pos[2] - state->z0), 2);
   vint += 0.5 * state->ktheta * pow((iTemp.theta - state->theta0), 2);
   eni += vint;
   return eni;
@@ -960,6 +1011,8 @@ void usage(){
                 "                   aLat = sqrt(3) hexSpace\n"
                 "                   bLat = hexSpace\n"
 		"   -k kz           Sets the value of kz in units of kb\n"
+		"   -t t           Sets the value of temperature in units of kelvin\n"
+		"   -q strength    Sets the strength of the dipole\n"
                 "\n",
                 
                 program_name);