anchor.diff:   
         Return
         END SUBROUTINE
  +
  
positive.code:       Subroutine Force(f,Upot)
        Implicit None
   
  !     Calculate The Forces And Potential Energy
         
        Include 'globals.inc'
  
        Integer:: I,J,Type1, Type2
        Double Precision:: Dx,Dy,Ff,R_square,R_square_i,R_six_i,Rcut, Rcutsq
        DOUBLE PRECISION::  K_bond, Rcut_bond,Rcutsq_bond 
  !**********Initialize the Forces, Potential Energy and Pressure to 0************
   
        DO I = 1,Natom
  
negative.code: 
  ! ROUGH DRAFT OF USING BONDS!
                  ! HEAD          TAIL1         TAIL2
              IF (Type1==2 .AND. Type2==3 .OR. Type2==4) THEN
               !A bond calculation required
  !QUESTION: should Energy_bond be Global?? Or, We can specify it here
  ! DONG: The Bonds I still don't quite understand but Energy_bond, bond_Rcut 
  ! should be global and defined from the read in of the data file. We 
  ! probably don't need to reassign them to K-bond, etc as I did here. 
  ! This is a rough draft. 
                K_bond = Energy_bond
                Rcut_bond = bond_Rcut
                Rcutsq_bond = Rcut_bond**2.0
                IF (R_square .Lt. Rcutsq_bond) THEN
                  sqR_square = sqrt(R_square)
                  Ff = -K_bond*(2.0-(Rcut_bond/sqR_square))
                  Upot = Upot + K_bond*((sqR_square - Rcut_bond)**2.0)  !The bond potential at Rcut = 0 
                  Press = Press + Ff
                  Fx(I) = Fx(I) + Ff*Dx
  