test_pot_bend.f90

! test_pot_bend.f90
! bend angle, cos(theta) potential
MODULE test_pot_module

  !------------------------------------------------------------------------------------------------!
  ! This software was written in 2016/17                                                           !
  ! by Michael P. Allen <m.p.allen@warwick.ac.uk>/<m.p.allen@bristol.ac.uk>                        !
  ! and Dominic J. Tildesley <d.tildesley7@gmail.com> ("the authors"),                             !
  ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"),     !
  ! published by Oxford University Press ("the publishers").                                       !
  !                                                                                                !
  ! LICENCE                                                                                        !
  ! Creative Commons CC0 Public Domain Dedication.                                                 !
  ! To the extent possible under law, the authors have dedicated all copyright and related         !
  ! and neighboring rights to this software to the PUBLIC domain worldwide.                        !
  ! This software is distributed without any warranty.                                             !
  ! You should have received a copy of the CC0 Public Domain Dedication along with this software.  !
  ! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.                               !
  !                                                                                                !
  ! DISCLAIMER                                                                                     !
  ! The authors and publishers make no warranties about the software, and disclaim liability       !
  ! for all uses of the software, to the fullest extent permitted by applicable law.               !
  ! The authors and publishers do not recommend use of this software for any purpose.              !
  ! It is made freely available, solely to clarify points made in the text. When using or citing   !
  ! the software, you should not imply endorsement by the authors or publishers.                   !
  !------------------------------------------------------------------------------------------------!

  USE, INTRINSIC :: iso_fortran_env, ONLY : error_unit

  IMPLICIT NONE
  PRIVATE

  ! Public routine
  PUBLIC :: force

  ! Public data
  INTEGER, PARAMETER, PUBLIC :: n = 3 ! Three-body potential

CONTAINS

  SUBROUTINE force  ( r, pot, f )
    IMPLICIT NONE
    REAL, DIMENSION(:,:),           INTENT(in)  :: r
    REAL,                           INTENT(out) :: pot
    REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: f

    ! We choose to make the polymer the minimum length needed for testing
    INTEGER                  :: a, b
    REAL, DIMENSION(3,2:n)   :: d     ! the d vectors
    REAL, DIMENSION(2:n,2:n) :: cc    ! the C coefficients
    REAL                     :: prefac, fac, fac1, fac2

    ! Routine to demonstrate the calculation of forces from the
    ! polymer angle-bending potential
    ! Written for ease of comparison with the text, rather than efficiency!

    ! Check dimensions to be sure
    IF ( ANY ( SHAPE(r) /= [3,n] ) ) THEN
       WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'r shape error', SHAPE(r), 3, n
       STOP 'Error in test_pot_bend'
    END IF

    ! Set up d vectors
    DO a = 2, n
       d(:,a) = r(:,a) - r(:,a-1)
    END DO

    ! Store C coefficients in a matrix
    ! In the general case we would not need to calculate every pair
    ! and also we would make use of the symmetry cc(b,a)=cc(a,b)
    DO a = 2, n
       DO b = 2, n
          cc(a,b) = DOT_PRODUCT ( d(:,a), d(:,b) )
       END DO
    END DO

    ! For this test there is just one angle
    a = n

    ! Here is the potential as a function of cos(theta)
    ! For testing we use the simplest form: v= -cos(theta)
    ! The notation matches that used in the appendix

    prefac = 1.0 /  SQRT(cc(a,a)*cc(a-1,a-1))
    fac    = cc(a,a-1)
    pot    = -prefac*fac ! This is -cos(theta)

    IF ( .NOT. PRESENT(f) ) RETURN

    IF ( ANY ( SHAPE(f) /= [3,n] ) ) THEN
       WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'f shape error', SHAPE(f), 3, n
       STOP 'Error in test_pot_bend'
    END IF

    ! Here we include the derivative of the potential with respect to cos(theta) in the prefactor
    ! For this simple case it is -1, so the forces are simply gradients of cos(theta) as in the text
    fac1 = fac / cc(a,a)
    fac2 = fac / cc(a-1,a-1)
    f(:,a)   = -prefac * ( fac1*d(:,a) - d(:,a-1) )
    f(:,a-1) =  prefac * ( fac1*d(:,a) - fac2*d(:,a-1) + d(:,a) - d(:,a-1) )
    f(:,a-2) =  prefac * ( fac2*d(:,a-1) - d(:,a) )

  END SUBROUTINE force

END MODULE test_pot_module