tortor_newmap Subroutine

public subroutine tortor_newmap(bds, d1, d2, ang1, ang2, v)

Uses

  • proc~~tortor_newmap~~UsesGraph proc~tortor_newmap tortor_newmap module~mod_utils mod_utils proc~tortor_newmap->module~mod_utils module~mod_memory mod_memory proc~tortor_newmap->module~mod_memory module~mod_utils->module~mod_memory module~mod_constants mod_constants module~mod_utils->module~mod_constants module~mod_memory->module~mod_constants iso_c_binding iso_c_binding module~mod_memory->iso_c_binding module~mod_io mod_io module~mod_memory->module~mod_io module~mod_constants->iso_c_binding module~mod_io->module~mod_constants

Store in module memory the data describing a new torsion-torsion map

Arguments

Type IntentOptional Attributes Name
type(ommp_bonded_type), intent(inout) :: bds
integer(kind=ip), intent(in) :: d1

Dimensions of the new map to be saved

integer(kind=ip), intent(in) :: d2

Dimensions of the new map to be saved

real(kind=rp), intent(in) :: ang1(:)

Value of torsion1 for the new map

real(kind=rp), intent(in) :: ang2(:)

Value of torsion2 for the new map

real(kind=rp), intent(in) :: v(:)

Value of potential for the new map


Calls

proc~~tortor_newmap~~CallsGraph proc~tortor_newmap tortor_newmap interface~mallocate mallocate proc~tortor_newmap->interface~mallocate interface~mfree mfree proc~tortor_newmap->interface~mfree proc~cyclic_spline cyclic_spline proc~tortor_newmap->proc~cyclic_spline proc~i_alloc1 i_alloc1 interface~mallocate->proc~i_alloc1 proc~i_alloc3 i_alloc3 interface~mallocate->proc~i_alloc3 proc~l_alloc2 l_alloc2 interface~mallocate->proc~l_alloc2 proc~r_alloc1 r_alloc1 interface~mallocate->proc~r_alloc1 proc~r_alloc3 r_alloc3 interface~mallocate->proc~r_alloc3 proc~i_alloc2 i_alloc2 interface~mallocate->proc~i_alloc2 proc~r_alloc2 r_alloc2 interface~mallocate->proc~r_alloc2 proc~l_alloc1 l_alloc1 interface~mallocate->proc~l_alloc1 proc~l_free1 l_free1 interface~mfree->proc~l_free1 proc~r_free3 r_free3 interface~mfree->proc~r_free3 proc~i_free2 i_free2 interface~mfree->proc~i_free2 proc~r_free2 r_free2 interface~mfree->proc~r_free2 proc~i_free1 i_free1 interface~mfree->proc~i_free1 proc~r_free1 r_free1 interface~mfree->proc~r_free1 proc~l_free2 l_free2 interface~mfree->proc~l_free2 proc~i_free3 i_free3 interface~mfree->proc~i_free3 proc~cyclic_spline->interface~mallocate proc~cyclic_spline->interface~mfree dgemm dgemm proc~cyclic_spline->dgemm dgetri dgetri proc~cyclic_spline->dgetri dgetrf dgetrf proc~cyclic_spline->dgetrf proc~chk_free chk_free proc~l_free1->proc~chk_free proc~r_free3->proc~chk_free proc~chk_alloc chk_alloc proc~i_alloc1->proc~chk_alloc proc~memory_init memory_init proc~i_alloc1->proc~memory_init proc~i_alloc3->proc~chk_alloc proc~i_alloc3->proc~memory_init proc~l_alloc2->proc~chk_alloc proc~l_alloc2->proc~memory_init proc~i_free2->proc~chk_free proc~r_free2->proc~chk_free proc~i_free1->proc~chk_free proc~r_alloc1->proc~chk_alloc proc~r_alloc1->proc~memory_init proc~r_alloc3->proc~chk_alloc proc~r_alloc3->proc~memory_init proc~i_alloc2->proc~chk_alloc proc~i_alloc2->proc~memory_init proc~r_alloc2->proc~chk_alloc proc~r_alloc2->proc~memory_init proc~l_alloc1->proc~chk_alloc proc~l_alloc1->proc~memory_init proc~r_free1->proc~chk_free proc~l_free2->proc~chk_free proc~i_free3->proc~chk_free proc~fatal_error fatal_error proc~chk_free->proc~fatal_error proc~chk_alloc->proc~fatal_error proc~ommp_message ommp_message proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output proc~close_output->proc~ommp_message

Called by

proc~~tortor_newmap~~CalledByGraph proc~tortor_newmap tortor_newmap proc~assign_tortors assign_tortors proc~assign_tortors->proc~tortor_newmap proc~mmpol_init_from_xyz mmpol_init_from_xyz proc~mmpol_init_from_xyz->proc~assign_tortors proc~ommp_system_from_qm_helper ommp_system_from_qm_helper proc~ommp_system_from_qm_helper->proc~assign_tortors proc~ommp_init_xyz ommp_init_xyz proc~ommp_init_xyz->proc~mmpol_init_from_xyz proc~c_ommp_system_from_qm_helper C_ommp_system_from_qm_helper proc~c_ommp_system_from_qm_helper->proc~ommp_system_from_qm_helper proc~c_ommp_init_xyz C_ommp_init_xyz proc~c_ommp_init_xyz->proc~ommp_init_xyz

Contents

Source Code


Source Code

    subroutine tortor_newmap(bds, d1, d2, ang1, ang2, v)
        !! Store in module memory the data describing a new torsion-torsion 
        !! map
        use mod_memory, only: mallocate, mfree
        use mod_utils, only: cyclic_spline

        implicit none

        type(ommp_bonded_type), intent(inout) :: bds
        ! Bonded potential data structure
        integer(ip), intent(in) :: d1, d2
        !! Dimensions of the new map to be saved
        real(rp), intent(in) :: ang1(:)
        !! Value of torsion1 for the new map 
        real(rp), intent(in) :: ang2(:)
        !! Value of torsion2 for the new map 
        real(rp), intent(in) :: v(:)
        !! Value of potential for the new map 

        integer :: i, j, ii
        real(rp), allocatable, dimension(:) :: a, b, c, d, dx, dy, dxy, tmpx, tmpy
        
        real(rp), allocatable :: rtmp(:)
        integer(ip), allocatable :: itmp(:,:)
        integer(ip) :: n_data, n_map

        if(allocated(bds%ttmap_ang1)) then
            ! Reallocate the arrays to make space for the new data
            n_data = size(bds%ttmap_ang1)
            call mallocate('torstors_newmap [rtmp]', n_data, rtmp)
            
            rtmp = bds%ttmap_ang1
            call mfree('torstors_newmap [ttmap_ang1]', bds%ttmap_ang1)
            call mallocate('torstors_newmap [ttmap_ang1]', &
                           n_data+d1*d2,  bds%ttmap_ang1)
            bds%ttmap_ang1(:n_data) = rtmp
            
            rtmp = bds%ttmap_ang2
            call mfree('torstors_newmap [ttmap_ang2]', bds%ttmap_ang2)
            call mallocate('torstors_newmap [ttmap_ang2]', &
                           n_data+d1*d2,  bds%ttmap_ang2)
            bds%ttmap_ang2(:n_data) = rtmp
            
            
            call mfree('torstors_newmap [rtmp]', rtmp)
            n_data = size(bds%ttmap_v)
            call mallocate('torstors_newmap [rtmp]', n_data, rtmp)
            
            rtmp = bds%ttmap_v
            call mfree('torstors_newmap [ttmap_v]', bds%ttmap_v)
            call mallocate('torstors_newmap [ttmap_v]', &
                           n_data+d1*d2,  bds%ttmap_v)
            bds%ttmap_v(:n_data) = rtmp
            
            rtmp = bds%ttmap_vx
            call mfree('torstors_newmap [ttmap_vx]', bds%ttmap_vx)
            call mallocate('torstors_newmap [ttmap_vx]', &
                           n_data+d1*d2,  bds%ttmap_vx)
            bds%ttmap_vx(:n_data) = rtmp

            rtmp = bds%ttmap_vy
            call mfree('torstors_newmap [ttmap_vy]', bds%ttmap_vy)
            call mallocate('torstors_newmap [ttmap_vy]', &
                           n_data+d1*d2,  bds%ttmap_vy)
            bds%ttmap_vy(:n_data) = rtmp

            rtmp = bds%ttmap_vxy
            call mfree('torstors_newmap [ttmap_vxy]', bds%ttmap_vxy)
            call mallocate('torstors_newmap [ttmap_vxy]', &
                           n_data+d1*d2,  bds%ttmap_vxy)
            bds%ttmap_vxy(:n_data) = rtmp
            call mfree('torstors_newmap [rtmp]', rtmp)

            n_map = size(bds%ttmap_shape, 2)
            call mallocate('torstors_newmap [itmp]', 2_ip, n_map, itmp)
            itmp = bds%ttmap_shape
            call mfree('torstors_newmap [ttmap_shape]', bds%ttmap_shape)
            call mallocate('torstors_newmap [ttmap_shape]', &
                           2_ip, n_map+1, bds%ttmap_shape)
            bds%ttmap_shape(:,:n_map) = itmp

            call mfree('torstors_newmap [itmp]', itmp)
        else 
            ! First allocation, n_data and n_map are just set for consistency
            n_data = 0
            n_map = 0
            call mallocate('torstors_newmap [ttmap_ang1]', d1*d2,  bds%ttmap_ang1)
            call mallocate('torstors_newmap [ttmap_ang2]', d1*d2,  bds%ttmap_ang2)
            call mallocate('torstors_newmap [ttmap_v]', d1*d2,  bds%ttmap_v)
            call mallocate('torstors_newmap [ttmap_vx]', d1*d2,  bds%ttmap_vx)
            call mallocate('torstors_newmap [ttmap_vy]', d1*d2,  bds%ttmap_vy)
            call mallocate('torstors_newmap [ttmap_vxy]', d1*d2,  bds%ttmap_vxy)
            call mallocate('torstors_newmap [ttmap_shape]', 2_ip, 1_ip, bds%ttmap_shape)
        end if

        call mallocate('tortor_newmap [a]', max(d1,d2), a)
        call mallocate('tortor_newmap [b]', max(d1,d2), b)
        call mallocate('tortor_newmap [c]', max(d1,d2), c)
        call mallocate('tortor_newmap [d]', max(d1,d2), d)
        call mallocate('tortor_newmap [dx]', d1*d2, dx)
        call mallocate('tortor_newmap [dy]', d1*d2, dy)
        call mallocate('tortor_newmap [dxy]', d1*d2, dxy)

        ! This part of the code computes df/dx, df/dy and d^2f/dxdy on the grid.
        ! Since we are basically interpolating on a sphere, we extract the 
        ! coordinate on a meridian, we interpolate it with a cubic spline, and
        ! finally we compute the derivative of this curve at the grid intersection
        ! The same is done in the second direction.
        ! To compute the mixed derivative we apply the same procedure but using
        ! the derivative data (basically we apply the procedure used to compute
        ! df/dx but using  df/dy data instead of actual f values.
        do i=1, d2
            call cyclic_spline(d1, ang1((i-1)*d1+1:i*d1), v((i-1)*d1+1:i*d1), &
                               a(1:d1), b(1:d1), c(1:d1), d(1:d1))
            dx((i-1)*d1+1:i*d1) = b(1:d1)
        end do
        
        ! df/dy since in this direction data are not contiguous, wa allocate 
        ! temporary arrays
        call mallocate('tortor_newmap [tmpx]', d2, tmpx)
        call mallocate('tortor_newmap [tmpy]', d2, tmpy)
        do i=1, d1
            ii = 1
            do j=i, (d2-1)*d1+i, d2
                tmpx(ii) = ang2(j)
                tmpy(ii) = v(j)
                ii = ii + 1
            end do
            call cyclic_spline(d2, tmpx, tmpy, &
                               a(1:d2), b(1:d2), c(1:d2), d(1:d2))
            
            ii = 1
            do j=i, (d2-1)*d1+i, d2
                dy(j) = b(ii)
                ii = ii + 1
            end do
        end do
        
        ! d^2f/dxdy in this case we use df/dx procedure to exploit data contiguity.
        do i=1, d2
            call cyclic_spline(d1, ang1((i-1)*d1+1:i*d1), dy((i-1)*d1+1:i*d1), &
                               a(1:d1), b(1:d1), c(1:d1), d(1:d1))
            dxy((i-1)*d1+1:i*d1) = b(1:d1)
        end do
        call mfree('tortor_newmap [tmpx]', tmpx)
        call mfree('tortor_newmap [tmpy]', tmpy)

        bds%ttmap_ang1(n_data+1:) = ang1
        bds%ttmap_ang2(n_data+1:) = ang2
        bds%ttmap_shape(1,n_map+1) = d1
        bds%ttmap_shape(2,n_map+1) = d2
        bds%ttmap_v(n_data+1:) = v
        bds%ttmap_vx(n_data+1:) = dx
        bds%ttmap_vy(n_data+1:) = dy
        bds%ttmap_vxy(n_data+1:) = dxy
        
        call mfree('tortor_newmap [a]', a)
        call mfree('tortor_newmap [b]', b)
        call mfree('tortor_newmap [c]', c)
        call mfree('tortor_newmap [d]', d)
        call mfree('tortor_newmap [dx]', dx)
        call mfree('tortor_newmap [dy]', dy)
        call mfree('tortor_newmap [dxy]', dxy)

    end subroutine tortor_newmap