cd $AMBERHOME
patch -p0 -N -r patch_rejects < NEB.patch 

Author: Daniel R. Roe
Date: 2010-03-25
Description:
Fixes several bugs in the NEB code. The neb_force array was not being zeroed
on end beads, which caused problems with certain compilers that do not 
automatically zero arrays. Also in certain places ncopy was being used instead
of neb_nbead, leading to incorrect energies being used in the tangent calc.
Some debugging code is added between #ifdef NEBDEBUG statements to make
future debugging easier.
Also updates the NEB test case runscripts to use up to 12 processors only,
test output from all beads, and use absolute error criterion instead of
truncating floating point digits.
-----------------------------------------------------------------------------
--- src/sander/pimd_force.f	2010-03-25 12:42:03.000000000 -0400
+++ src/sander/pimd_force.f	2010-03-25 12:23:14.000000000 -0400
@@ -457,17 +457,21 @@
 !  i.e. first and last still have both neighbors defined.
 !**********************************************************
 
-  ! these ranks imply we use commmaster for the communicator
-  next_node = masterrank+1
-  if(next_node>=mastersize) next_node = 0
-  prev_node = masterrank-1
-  if(prev_node<0) prev_node = mastersize-1
+   ! these ranks imply we use commmaster for the communicator
+   next_node = masterrank+1
+   if(next_node>=mastersize) next_node = 0
+   prev_node = masterrank-1
+   if(prev_node<0) prev_node = mastersize-1
  
+#ifdef NEBDEBUG
+   write(50+worldrank,'(a,i3,a,i3)') "NEBDBG: next_node=",next_node," prev_node=",prev_node
+   write(50+worldrank,'(a,i8)') "NEBDBG: last_neb_atom=",last_neb_atom
+   write(50+worldrank,'(a,i8)') "NEBDBG: nattgtrms=",nattgtrms
+   write(50+worldrank,'(a,i3)') "NEBDBG: neb_nbead=",neb_nbead
+#endif
 
-
-
-      ! carlos: only send/receive neighbor coordinates!
-      ! carlos: beads 1 and neb_nbead should just exit! (no coordinate update, and no neighbor)
+   ! carlos: only send/receive neighbor coordinates!
+   ! carlos: beads 1 and neb_nbead should just exit! (no coordinate update, and no neighbor)
 
       if (mod(mybeadid,2)==0) then  !even, exchange with bead+1
 
@@ -499,6 +503,10 @@
 
       endif
 
+#ifdef NEBDEBUG
+   write(50+worldrank,'(a,f12.4)') "NEBDBG: epot=",epot
+#endif
+
       call mpi_allgather(epot,1,MPI_DOUBLE_PRECISION, &
                          neb_nrg_all,1, MPI_DOUBLE_PRECISION, &
                          commmaster,ierr)
@@ -514,13 +522,26 @@
         call rmsfit( xnext, x, mass, fitgp,rmsgp, rmsdvalue, nattgtrms, nattgtfit, rmsok) 
 
    !ezero is the higher energy of the two end points
-   ezero = max( neb_nrg_all(1), neb_nrg_all(ncopy) )
+#ifdef NEBDEBUG
+   write(50+worldrank,'(a,i3)') "NEBDBG: Getting max (ezero) between 1 and ",neb_nbead
+#endif
+   ezero = max( neb_nrg_all(1), neb_nrg_all(neb_nbead) )
+#ifdef NEBDEBUG
+   write(50+worldrank,'(a,f12.4)') "NEBDBG: ezero=",ezero
+   write(50+worldrank,'(a,f12.4)') "NEBDBG:   neb_nrg_all(1)=",neb_nrg_all(1)
+#endif
    
    !emax is the highest energy of all the points
    emax = neb_nrg_all(1)
-   do rep = 2, ncopy
+   do rep = 2, neb_nbead
+#ifdef NEBDEBUG
+      write(50+worldrank,'(a,i3,a,f12.4)') "NEBDBG:   neb_nrg_all(",rep,")=",neb_nrg_all(rep)
+#endif
       emax = max(emax, neb_nrg_all(rep))
    end do
+#ifdef NEBDEBUG
+   write(50+worldrank,'(a,f12.4)') "NEBDBG: emax=",emax
+#endif
    
    !spring2 is the spring constant of the second point in path
    if (skmin.EQ.skmax) then
@@ -539,13 +560,24 @@
    !nebrms = 0.d0
 
 ! don't bother with forces for first or last bead (forces are zeroed anyway)
+! DAN ROE: Zero out the forces on all beads here since first and last beads
+!          are about to exit. Some compilers do NOT automatically zero
+!          arrays. Call in force.f only uses nattgtrms*3 spots
+   do i=1,nattgtrms*3
+     neb_force(i)=0.d0
+   enddo
 
 ! CHANGE LATER IF WE DON'T DO (FAKE) MD ON ENDPOINTS!
    if( mybeadid==1.or.mybeadid==neb_nbead) then
         !write (6,*) "I am first or last bead, returning "
+#ifdef NEBDEBUG
+      write(worldrank+50,'(a)') "-------------------------------------------------"
+#endif
         return
    endif
-
+#ifdef NEBDEBUG
+   write(50+worldrank,'(a,i3,a,i3,a,i3)') "NEBDBG: Rep=",rep," rep+1=",rep+1," rep-1=",rep-1
+#endif
    !calculate spring constant for rep and rep+1
  
    spring = spring2
@@ -558,7 +590,9 @@
       spring2 = skmax - skmin
    end if
 
-   tangents = 0.0  
+   !tangents = 0.0 
+! DAN NOTE: tangents is set to 0.0 here, also in the following do loop. Is this
+! kind of assignment ok? 
 
 ! carlos: how are rep+1 for last rep and rep-1 for first handled as indices??
 ! for now it doesn't matter since first and last rep exited already
@@ -636,6 +670,15 @@
       end do
    end if
 
+#ifdef NEBDEBUG
+   do i=1,nattgtrms
+      iatom=rmsgp(i)
+      iatm3 = 3*(iatom - 1)
+      write(worldrank+50,'(a,f12.4,a,f12.4,a,f12.4)') &
+         "NEBDBG: t1=",tangents(iatm3+1)," t2=",tangents(iatm3+2)," t3=",tangents(iatm3+3)
+   enddo
+#endif
+
    call normalize(tangents, 3 * last_neb_atom)
 
 
@@ -725,6 +768,10 @@
 
    !nebrms = dotproduct
 
+#ifdef NEBDEBUG
+   write(worldrank+50,'(a)') "-------------------------------------------------"
+#endif
+
     return
 end subroutine full_neb_forces
 
--- test/neb-testcases/neb_explicit/Run.neb_explicit	2010-03-25 12:44:27.000000000 -0400
+++ test/neb-testcases/neb_explicit/Run.neb_explicit	2010-03-25 12:47:57.000000000 -0400
@@ -13,17 +13,17 @@
   exit(0)
 else
   set numprocs=`echo $DO_PARALLEL | awk -f ../../numprocs.awk `
-  if ( $numprocs == 4 || $numprocs == 8 || $numprocs == 16 || $numprocs == 24 ) then
+  if ( $numprocs == 4 || $numprocs == 8 || $numprocs == 12 ) then
       goto runtest
   else if ( $?MP_PROCS) then
-      if ( $MP_PROCS == 4 || $MP_PROCS == 8 || $MP_PROCS == 16 || $MP_PROCS == 24) then
+      if ( $MP_PROCS == 4 || $MP_PROCS == 8 || $MP_PROCS == 12 ) then
         goto runtest
       endif
   endif
 endif
 
 echo " This test case requires a least 4 mpi threads."
-echo " The number of mpi threads must also be a multiple of 4 and not more than 24."
+echo " The number of mpi threads must also be a multiple of 4 and not more than 12."
 echo " Not running test, exiting....."
 exit(0)
 
@@ -52,7 +52,10 @@
 EOF
 
 $DO_PARALLEL $sander -ng 4 -groupfile groupfile.in || goto error 
-../../dacdif -t 1 neb_explicit_01.out.save neb_explicit_01.out
+../../dacdif neb_explicit_01.out.save neb_explicit_01.out -a 1E-3
+../../dacdif neb_explicit_02.out.save neb_explicit_02.out -a 1E-3
+../../dacdif neb_explicit_03.out.save neb_explicit_03.out -a 1E-3
+../../dacdif neb_explicit_04.out.save neb_explicit_04.out -a 1E-3
 
 /bin/rm -f mdin *.inf *.mdcrd *.rst *.out 
 endif
--- test/neb-testcases/neb_gb_full/Run.neb_gb_full	2010-03-25 12:45:12.000000000 -0400
+++ test/neb-testcases/neb_gb_full/Run.neb_gb_full	2010-03-25 12:41:03.000000000 -0400
@@ -13,17 +13,17 @@
   exit(0)
 else
   set numprocs=`echo $DO_PARALLEL | awk -f ../../numprocs.awk `
-  if ( $numprocs == 4 || $numprocs == 8 || $numprocs == 16 || $numprocs == 24 ) then
+  if ( $numprocs == 4 || $numprocs == 8 || $numprocs == 12 ) then
       goto runtest
   else if ( $?MP_PROCS)then
-      if ( $MP_PROCS == 4 || $MP_PROCS == 8 || $MP_PROCS == 16 || $MP_PROCS == 24)then
+      if ( $MP_PROCS == 4 || $MP_PROCS == 8 || $MP_PROCS == 12 )then
         goto runtest
       endif
   endif
 endif
 
 echo " This test case requires a least 4 mpi threads."
-echo " The number of mpi threads must also be a multiple of 4 and not more than 24."
+echo " The number of mpi threads must also be a multiple of 4 and not more than 12."
 echo " Not running test, exiting....."
 exit(0)
 
@@ -53,7 +53,10 @@
 EOF
 
 $DO_PARALLEL $sander -ng 4 -groupfile groupfile.in || goto error
-../../dacdif -t 1 neb_gb_full_01.out.save neb_gb_full_01.out
+../../dacdif neb_gb_full_01.out.save neb_gb_full_01.out -a 1E-3
+../../dacdif neb_gb_full_02.out.save neb_gb_full_02.out -a 1E-3
+../../dacdif neb_gb_full_03.out.save neb_gb_full_03.out -a 1E-3
+../../dacdif neb_gb_full_04.out.save neb_gb_full_04.out -a 1E-3
 
 /bin/rm -f mdin *.inf *.mdcrd *.rst *.out 
 endif
--- test/neb-testcases/neb_gb_partial/Run.neb_gb_partial	2010-03-25 12:44:47.000000000 -0400
+++ test/neb-testcases/neb_gb_partial/Run.neb_gb_partial	2010-03-25 12:40:31.000000000 -0400
@@ -13,17 +13,17 @@
   exit(0)
 else
   set numprocs=`echo $DO_PARALLEL | awk -f ../../numprocs.awk `
-  if ( $numprocs == 4 || $numprocs == 8 || $numprocs == 16 || $numprocs == 24 ) then
+  if ( $numprocs == 4 || $numprocs == 8 || $numprocs == 12 ) then
       goto runtest
   else if ( $?MP_PROCS)then
-      if ( $MP_PROCS == 4 || $MP_PROCS == 8 || $MP_PROCS == 16 || $MP_PROCS == 24)then
+      if ( $MP_PROCS == 4 || $MP_PROCS == 8 || $MP_PROCS == 12 )then
         goto runtest
       endif
   endif
 endif
 
 echo " This test case requires a least 4 mpi threads."
-echo " The number of mpi threads must also be a multiple of 4 and not more than 24."
+echo " The number of mpi threads must also be a multiple of 4 and not more than 12."
 echo " Not running test, exiting....."
 exit(0)
 
@@ -53,7 +53,10 @@
 EOF
 
 $DO_PARALLEL $sander -ng 4 -groupfile groupfile.in || goto error
-../../dacdif -t 1 neb_gb_partial_01.out.save neb_gb_partial_01.out
+../../dacdif neb_gb_partial_01.out.save neb_gb_partial_01.out -a 1E-3
+../../dacdif neb_gb_partial_02.out.save neb_gb_partial_02.out -a 1E-3
+../../dacdif neb_gb_partial_03.out.save neb_gb_partial_03.out -a 1E-3
+../../dacdif neb_gb_partial_04.out.save neb_gb_partial_04.out -a 1E-3
 
 /bin/rm -f mdin *.inf *.mdcrd *.rst *.out 
 endif