/src/harmonicrestraintenergy.F

http://protoms.googlecode.com/ · FORTRAN Legacy · 79 lines · 52 code · 10 blank · 17 comment · 6 complexity · d2f40fda9b5ac37517602062c148d55a MD5 · raw file 

    

  1. double precision function harmonicRestraintEnergy(savenrg,inrgb,inrgf)

    2. 

  2.       include 'dimensions.inc'

    3. 

  3.       include 'movelist.inc'

    4. 

  4.       include 'printer.inc'

    5. 

  5.       include 'constants.inc'

    6. 

  6.       include 'enums.inc'

    7. 

  7.       include 'extraenergies.inc'

    8. 

  8. c###############################################################################

    9. 

  9. c

    10. 

  10. c     This function calculates the total harmonic restraint energy

    11. 

  11. c     for all of the atoms that have moved.

    12. 

  12. c

    13. 

  13. c     (C) Christopher Woods, November 2004

    14. 

  14. c

    15. 

  15. c###############################################################################

    16. 

  16. 

    17. 

  17.       integer savenrg

    18. 

  18.       integer i

    19. 

  19.       double precision inrgb,inrgf

    20. 

  20.       double precision ic(3),icb(3),icf(3),dist,distb,distf

    21. 

  21.       logical hasAtomMoved

    22. 

  22.       

    23. 

  23.       harmonicRestraintEnergy = ZERO

    24. 

  24.       inrgb = ZERO

    25. 

  25.       inrgf = ZERO

    26. 

  26.       

    27. 

  27.       if (nHarmonicRestraints.le.0) return

    28. 

  28.       

    29. 

  29.       if (allmoved) then

    30. 

  30. c       everything has moved - evaluate all of the harmonic energies

    31. 

  31.         do i=1,nHarmonicRestraints

    32. 

  32. c         get the coordinates of the restrained atom

    33. 

  33.           call getCoordinates(HarmAtom(i,1),HarmAtom(i,2),HarmAtom(i,3),ic,icb,icf)

    34. 

  34. c         calculate the distance to the restrained point

    35. 

  35.           dist = (ic(1)-HarmPoint(i,1))**2 + (ic(2)-HarmPoint(i,2))**2 + (ic(3)-HarmPoint(i,3))**2

    36. 

  36.           distb = (icb(1)-HarmPoint(i,1))**2 + (icb(2)-HarmPoint(i,2))**2 + (icb(3)-HarmPoint(i,3))**2

    37. 

  37.           distf = (icf(1)-HarmPoint(i,1))**2 + (icf(2)-HarmPoint(i,2))**2 + (icf(3)-HarmPoint(i,3))**2

    38. 

  38.           

    39. 

  39.           harmonicRestraintEnergy = harmonicRestraintEnergy + HarmConst(i)*dist

    40. 

  40.           inrgb = inrgb + HarmConst(i)*distb

    41. 

  41.           inrgf = inrgf + HarmConst(i)*distf

    42. 

  42.         enddo

    43. 

  43.       else

    44. 

  44. c       loop through all of the restraints and see if they need updating

    45. 

  45.         do i=1,nHarmonicRestraints

    46. 

  46.           if (hasAtomMoved(HarmAtom(i,1),HarmAtom(i,2),HarmAtom(i,3))) then

    47. 

  47.             call getCoordinates(HarmAtom(i,1),HarmAtom(i,2),HarmAtom(i,3),ic,icb,icf)

    48. 

  48. 

    49. 

  49. c           calculate the distance to the restrained point

    50. 

  50.             dist = (ic(1)-HarmPoint(i,1))**2 + (ic(2)-HarmPoint(i,2))**2 + (ic(3)-HarmPoint(i,3))**2

    51. 

  51.             distb = (icb(1)-HarmPoint(i,1))**2 + (icb(2)-HarmPoint(i,2))**2 + (icb(3)-HarmPoint(i,3))**2

    52. 

  52.             distf = (icf(1)-HarmPoint(i,1))**2 + (icf(2)-HarmPoint(i,2))**2 + (icf(3)-HarmPoint(i,3))**2

    53. 

  53.           

    54. 

  54.             harmonicRestraintEnergy = harmonicRestraintEnergy + HarmConst(i)*dist

    55. 

  55.             inrgb = inrgb + HarmConst(i)*distb

    56. 

  56.             inrgf = inrgf + HarmConst(i)*distf

    57. 

  57.           endif

    58. 

  58.         enddo

    59. 

  59.       endif      

    60. 

  60.         

    61. 

  61.       if (savenrg.eq.OLD) then

    62. 

  62.         oldHarmEnergyPart = harmonicRestraintEnergy

    63. 

  63.         oldHarmEnergyPartB = inrgb

    64. 

  64.         oldHarmEnergyPartF = inrgf

    65. 

  65.         write(printstring,*) 'old harmNrg ?',oldHarmEnergyPart,' F ',oldHarmEnergyPartF,

    66. 

  66.      .    ' B ',oldHarmEnergyPartB

    67. 

  67.         call printLine(DEBUG,printstring)      

    68. 

  68.       else if (savenrg.eq.NEW) then

    69. 

  69.         newHarmEnergyPart = harmonicRestraintEnergy

    70. 

  70.         newHarmEnergyPartB = inrgb

    71. 

  71.         newHarmEnergyPartF = inrgf

    72. 

  72.         write(printstring,*) 'new harmNrg ?',newHarmEnergyPart,' F ',newHarmEnergyPartF,

    73. 

  73.      .    ' B ',newHarmEnergyPartB

    74. 

  74.         call printLine(DEBUG,printstring)      

    75. 

  75.       endif      

    76. 

  76. 

    77. 

  77.       return

    78. 

  78.       end
    79.