angtor_potential Subroutine

public subroutine angtor_potential(bds, V)

Arguments

Type IntentOptional Attributes Name
type(ommp_bonded_type), intent(in) :: bds
real(kind=rp), intent(inout) :: V

Called by

proc~~angtor_potential~~CalledByGraph proc~angtor_potential angtor_potential proc~ommp_get_angtor_energy ommp_get_angtor_energy proc~ommp_get_angtor_energy->proc~angtor_potential proc~ommp_get_full_bnd_energy ommp_get_full_bnd_energy proc~ommp_get_full_bnd_energy->proc~angtor_potential proc~c_ommp_get_angtor_energy C_ommp_get_angtor_energy proc~c_ommp_get_angtor_energy->proc~ommp_get_angtor_energy program~test_si_potential test_SI_potential program~test_si_potential->proc~ommp_get_angtor_energy proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_bnd_energy proc~c_ommp_get_full_bnd_energy C_ommp_get_full_bnd_energy proc~c_ommp_get_full_bnd_energy->proc~ommp_get_full_bnd_energy proc~c_ommp_get_full_energy C_ommp_get_full_energy proc~c_ommp_get_full_energy->proc~ommp_get_full_energy

Contents

Source Code


Source Code

    subroutine angtor_potential(bds, V)

        implicit none

        type(ommp_bonded_type), intent(in) :: bds
        ! Bonded potential data structure
        real(rp), intent(inout) :: V
        real(rp) :: thet, costhet, dihef(3), delta_a(2), vat, l1, l2, &
                    dr1(3), dr2(3), angle1, angle2
        integer(ip) :: i, j, k, ia1, ia2
        
        if(.not. bds%use_angtor) return

        !$omp parallel do default(shared) reduction(+:V) &
        !$omp private(i,costhet,thet,j,dihef,ia1,ia2,dr1,dr2,l1,l2,angle1,angle2,delta_a,vat,k)
        do i=1, bds%nangtor
            ! Atoms that defines the dihedral angle
            costhet = cos_torsion(bds%top, bds%angtorat(:,i))
            thet = acos(costhet)
            do j=1, 3
                dihef(j) = 1.0 + cos(j*thet+bds%torsphase(j,bds%angtor_t(i)))
            end do

            ia1 = bds%angtor_a(1,i)
            ia2 = bds%angtor_a(2,i)
            
            dr1 = bds%top%cmm(:, bds%angleat(1,ia1)) - &
                  bds%top%cmm(:, bds%angleat(2,ia1))
            dr2 = bds%top%cmm(:, bds%angleat(3,ia1)) - &
                  bds%top%cmm(:, bds%angleat(2,ia1))
            l1 = norm2(dr1)
            l2 = norm2(dr2)
            angle1 = acos(dot_product(dr1, dr2)/(l1*l2))

            dr1 = bds%top%cmm(:, bds%angleat(1,ia2)) - &
                  bds%top%cmm(:, bds%angleat(2,ia2))
            dr2 = bds%top%cmm(:, bds%angleat(3,ia2)) - &
                  bds%top%cmm(:, bds%angleat(2,ia2))
            l1 = norm2(dr1)
            l2 = norm2(dr2)
            angle2 = acos(dot_product(dr1, dr2)/(l1*l2))
           
            delta_a(1) = angle1 - bds%eqangle(bds%angtor_a(1,i))
            delta_a(2) = angle2 - bds%eqangle(bds%angtor_a(2,i))

            do j=1,2
                vat = 0.0
                do k=1, 3
                    vat = vat + bds%angtork((j-1)*3+k,i) * dihef(k)
                end do
                V = V + vat * delta_a(j)
            end do
        end do

    end subroutine angtor_potential