assign_vdw Subroutine

public subroutine assign_vdw(vdw, top, prm_buf)

Uses

  • proc~~assign_vdw~~UsesGraph proc~assign_vdw assign_vdw module~mod_nonbonded mod_nonbonded proc~assign_vdw->module~mod_nonbonded module~mod_constants mod_constants proc~assign_vdw->module~mod_constants module~mod_io mod_io proc~assign_vdw->module~mod_io module~mod_memory mod_memory proc~assign_vdw->module~mod_memory module~mod_nonbonded->module~mod_constants module~mod_nonbonded->module~mod_memory module~mod_adjacency_mat mod_adjacency_mat module~mod_nonbonded->module~mod_adjacency_mat module~mod_topology mod_topology module~mod_nonbonded->module~mod_topology module~mod_neighbor_list mod_neighbor_list module~mod_nonbonded->module~mod_neighbor_list iso_c_binding iso_c_binding module~mod_constants->iso_c_binding module~mod_io->module~mod_constants module~mod_memory->module~mod_constants module~mod_memory->module~mod_io module~mod_memory->iso_c_binding module~mod_adjacency_mat->module~mod_memory module~mod_topology->module~mod_memory module~mod_topology->module~mod_adjacency_mat module~mod_neighbor_list->module~mod_io module~mod_neighbor_list->module~mod_memory module~mod_neighbor_list->module~mod_adjacency_mat

Arguments

Type IntentOptional Attributes Name
type(ommp_nonbonded_type), intent(inout) :: vdw

Non-bonded structure to be initialized
type(ommp_topology_type), intent(inout) :: top

Topology structure
character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)

Char buffer containing the prm file loaded in RAM

Calls

proc~~assign_vdw~~CallsGraph proc~assign_vdw assign_vdw proc~read_atom_cards read_atom_cards proc~assign_vdw->proc~read_atom_cards interface~mallocate mallocate proc~assign_vdw->interface~mallocate proc~tokenize tokenize proc~assign_vdw->proc~tokenize proc~isreal isreal proc~assign_vdw->proc~isreal proc~fatal_error fatal_error proc~assign_vdw->proc~fatal_error proc~isint isint proc~assign_vdw->proc~isint proc~vdw_init vdw_init proc~assign_vdw->proc~vdw_init proc~vdw_set_pair vdw_set_pair proc~assign_vdw->proc~vdw_set_pair interface~mfree mfree proc~assign_vdw->interface~mfree proc~read_atom_cards->interface~mallocate proc~read_atom_cards->proc~tokenize proc~read_atom_cards->proc~fatal_error proc~read_atom_cards->proc~isint proc~read_atom_cards->interface~mfree proc~count_substr_occurence count_substr_occurence proc~read_atom_cards->proc~count_substr_occurence 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_alloc3 r_alloc3 interface~mallocate->proc~r_alloc3 proc~i_alloc2 i_alloc2 interface~mallocate->proc~i_alloc2 proc~r_alloc1 r_alloc1 interface~mallocate->proc~r_alloc1 proc~r_alloc2 r_alloc2 interface~mallocate->proc~r_alloc2 proc~l_alloc1 l_alloc1 interface~mallocate->proc~l_alloc1 proc~ommp_message ommp_message proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output proc~vdw_init->interface~mallocate proc~vdw_init->proc~fatal_error proc~nl_init nl_init proc~vdw_init->proc~nl_init proc~vdw_set_pair->interface~mallocate proc~vdw_set_pair->interface~mfree proc~l_free1 l_free1 interface~mfree->proc~l_free1 proc~i_free2 i_free2 interface~mfree->proc~i_free2 proc~r_free3 r_free3 interface~mfree->proc~r_free3 proc~r_free2 r_free2 interface~mfree->proc~r_free2 proc~l_free2 l_free2 interface~mfree->proc~l_free2 proc~i_free1 i_free1 interface~mfree->proc~i_free1 proc~r_free1 r_free1 interface~mfree->proc~r_free1 proc~i_free3 i_free3 interface~mfree->proc~i_free3 proc~chk_free chk_free proc~l_free1->proc~chk_free proc~i_free2->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~r_free2->proc~chk_free proc~close_output->proc~ommp_message proc~l_free2->proc~chk_free proc~i_free1->proc~chk_free proc~r_free1->proc~chk_free 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_alloc1->proc~chk_alloc proc~r_alloc1->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~nl_init->interface~mallocate proc~nl_init->proc~fatal_error proc~nl_update nl_update proc~nl_init->proc~nl_update proc~i_free3->proc~chk_free proc~chk_free->proc~fatal_error proc~chk_alloc->proc~fatal_error proc~nl_update->interface~mfree proc~nl_update->proc~ommp_message proc~time_push time_push proc~nl_update->proc~time_push proc~time_pull time_pull proc~nl_update->proc~time_pull proc~reverse_grp_tab reverse_grp_tab proc~nl_update->proc~reverse_grp_tab proc~time_push->proc~fatal_error proc~mem_stat mem_stat proc~time_push->proc~mem_stat proc~time_pull->proc~fatal_error proc~time_pull->proc~ommp_message proc~time_pull->proc~mem_stat proc~reverse_grp_tab->interface~mallocate proc~reverse_grp_tab->interface~mfree proc~compress_list compress_list proc~reverse_grp_tab->proc~compress_list proc~mem_stat->proc~memory_init proc~compress_list->interface~mallocate proc~compress_list->interface~mfree

Called by

proc~~assign_vdw~~CalledByGraph proc~assign_vdw assign_vdw proc~mmpol_init_from_xyz mmpol_init_from_xyz proc~mmpol_init_from_xyz->proc~assign_vdw proc~qm_helper_init_vdw_prm qm_helper_init_vdw_prm proc~qm_helper_init_vdw_prm->proc~assign_vdw proc~ommp_system_from_qm_helper ommp_system_from_qm_helper proc~ommp_system_from_qm_helper->proc~assign_vdw 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_qm_helper_init_vdw_prm C_ommp_qm_helper_init_vdw_prm proc~c_ommp_qm_helper_init_vdw_prm->proc~qm_helper_init_vdw_prm proc~ommp_create_link_atom ommp_create_link_atom proc~ommp_create_link_atom->proc~qm_helper_init_vdw_prm proc~c_ommp_init_xyz C_ommp_init_xyz proc~c_ommp_init_xyz->proc~ommp_init_xyz proc~c_ommp_create_link_atom C_ommp_create_link_atom proc~c_ommp_create_link_atom->proc~ommp_create_link_atom

Contents

Source Code

assign_vdw

Source Code

    subroutine assign_vdw(vdw, top, prm_buf)
        use mod_memory, only: mallocate, mfree
        use mod_io, only: fatal_error
        use mod_nonbonded, only: ommp_nonbonded_type, vdw_init, vdw_set_pair
        use mod_constants, only: angstrom2au, kcalmol2au, OMMP_DEFAULT_NL_CUTOFF
        
        implicit none
       
        type(ommp_nonbonded_type), intent(inout) :: vdw
        !! Non-bonded structure to be initialized
        type(ommp_topology_type), intent(inout) :: top
        !! Topology structure
        character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)
        !! Char buffer containing the prm file loaded in RAM

        integer(ip) :: il, i, j, l, tokb, toke
        character(len=OMMP_STR_CHAR_MAX) :: line, errstring
        character(len=20) :: radrule, radsize, radtype, vdwtype, epsrule
        integer(ip), allocatable :: vdwat(:), vdwpr_a(:), vdwpr_b(:)
        real(rp), allocatable :: vdw_e_prm(:), vdw_r_prm(:), vdw_f_prm(:), &
                                 vdwpr_r(:), vdwpr_e(:)
        integer(ip) :: nvdw, ivdw, atc, nvdwpr, ivdwpr
        logical :: done
        logical(lp), allocatable :: maska(:), maskb(:)

        if(.not. top%atclass_initialized .or. .not. top%atz_initialized) then
            call read_atom_cards(top, prm_buf)
        end if
        
        ! Read all the lines of file just to count how large vector should be 
        ! allocated 
        nvdw = 0
        nvdwpr = 0
        do il=1, size(prm_buf) 
            line = prm_buf(il)
            if(line(:4) == 'vdw ') nvdw = nvdw + 1
            if(line(:6) == 'vdwpr ' .or. line(:8) == 'vdwpair ') &
                nvdwpr = nvdwpr + 1
        end do

        ! VDW
        call mallocate('read_prm [vdwat]', nvdw, vdwat)
        call mallocate('read_prm [vdw_r_prm]', nvdw, vdw_r_prm)
        call mallocate('read_prm [vdw_e_prm]', nvdw, vdw_e_prm)
        call mallocate('read_prm [vdw_f_prm]', nvdw, vdw_f_prm)
        vdw_f_prm = 1.0_rp
        ivdw = 1

        ! VDWPR
        call mallocate('read_prm [vdwpr_a]', nvdwpr, vdwpr_a)
        call mallocate('read_prm [vdwpr_b]', nvdwpr, vdwpr_b)
        call mallocate('read_prm [vdwpr_e]', nvdwpr, vdwpr_e)
        call mallocate('read_prm [vdwpr_r]', nvdwpr, vdwpr_r)
        allocate(maska(top%mm_atoms), maskb(top%mm_atoms))
        ivdwpr = 1

        ! Default rules for VDW (from Tinker Manual)
        radrule = "arithmetic"
        radsize = "radius"
        radtype = "r-min"
        vdwtype = "lennard-jones"
        epsrule = "geometric"
        ! Those are defaults defined in Tinker manual
        vdw%vdw_screening = [0.0, 0.0, 1.0, 1.0]

        ! Restart the reading from the beginning to actually save the parameters
        i=1
        do il=1, size(prm_buf) 
            line = prm_buf(il)
           
            if(line(:13) == 'vdw-12-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW-12-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdw%vdw_screening(1)
                if(vdw%vdw_screening(1) > 1.0) &
                    vdw%vdw_screening(1) = 1/vdw%vdw_screening(1)
            
            else if(line(:13) == 'vdw-13-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW-13-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdw%vdw_screening(2)
                if(vdw%vdw_screening(2) > 1.0) &
                    vdw%vdw_screening(2) = 1/vdw%vdw_screening(2)
            
            else if(line(:13) == 'vdw-14-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW-14-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdw%vdw_screening(3)
                if(vdw%vdw_screening(3) > 1.0) &
                    vdw%vdw_screening(3) = 1/vdw%vdw_screening(3)
            
            else if(line(:13) == 'vdw-15-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW-15-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdw%vdw_screening(4)
                if(vdw%vdw_screening(4) > 1.0) &
                    vdw%vdw_screening(4) = 1/vdw%vdw_screening(4)

            else if(line(:12) == 'epsilonrule ') then
                tokb = 12
                toke = tokenize(line, tokb)
                read(line(tokb:toke), *) epsrule
            
            else if(line(:8) == 'vdwtype ') then
                tokb = 9
                toke = tokenize(line, tokb)
                read(line(tokb:toke), *) vdwtype
            
            else if(line(:11) == 'radiusrule ') then
                tokb = 12
                toke = tokenize(line, tokb)
                read(line(tokb:toke), *) radrule

            else if(line(:11) == 'radiussize ') then
                tokb = 12
                toke = tokenize(line, tokb)
                read(line(tokb:toke), *) radsize

            else if(line(:11) == 'radiustype ') then
                tokb = 12
                toke = tokenize(line, tokb)
                read(line(tokb:toke), *) radtype

            else if(line(:4) == 'vdw ') then
                tokb = 5
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdwat(ivdw)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdw_r_prm(ivdw)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDW card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdw_e_prm(ivdw)
                
                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(toke > 0) then
                    if(.not. isreal(line(tokb:toke))) then
                        write(errstring, *) "Wrong VDW card"
                        call fatal_error(errstring)
                    end if
                    read(line(tokb:toke), *) vdw_f_prm(ivdw)
                endif

                ivdw = ivdw + 1
            else if(line(:6) == 'vdwpr ' .or. line(:8) == 'vdwpair ') then
                if(line(:6) == 'vdwpr ') then
                    tokb = 7
                else
                    tokb = 9
                end if
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDWPR card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdwpr_a(ivdwpr)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDWPR card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdwpr_b(ivdwpr)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDWPR card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdwpr_r(ivdwpr)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong VDWPR card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) vdwpr_e(ivdwpr)

                ivdwpr = ivdwpr + 1
            end if
            i = i+1
        end do
        
        call vdw_init(vdw, top, vdwtype, radrule, radsize, &
                      radtype, epsrule, OMMP_DEFAULT_NL_CUTOFF)
        
        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,atc,done)
        do i=1, top%mm_atoms
            ! Atom class for current atom
            atc = top%atclass(i)
            
            ! VdW parameters
            done = .false.
            do j=1, nvdw
                if(vdwat(j) == atc) then
                    done = .true.
                    vdw%vdw_e(i) = vdw_e_prm(j) * kcalmol2au
                    vdw%vdw_r(i) = vdw_r_prm(j) * angstrom2au
                    vdw%vdw_f(i) = vdw_f_prm(j)
                    exit
                end if
            end do
            if(.not. done) then
                call fatal_error("VdW parameter not found!")
            end if
        end do
            
        ! VdW pair parameters
        do l=1, nvdwpr
            maska = (top%atclass == vdwpr_a(l))
            maskb = (top%atclass == vdwpr_b(l))
            call vdw_set_pair(vdw, maska, maskb, &
                              vdwpr_r(l) * angstrom2au, &
                              vdwpr_e(l) * kcalmol2au)
        end do
        
        call mfree('read_prm [vdwat]', vdwat)
        call mfree('read_prm [vdw_r_prm]', vdw_r_prm)
        call mfree('read_prm [vdw_e_prm]', vdw_e_prm)
        call mfree('read_prm [vdw_f_prm]', vdw_f_prm)
        call mfree('read_prm [vdwpr_a]', vdwpr_a)
        call mfree('read_prm [vdwpr_b]', vdwpr_b)
        call mfree('read_prm [vdwpr_e]', vdwpr_e)
        call mfree('read_prm [vdwpr_r]', vdwpr_r)
        deallocate(maska, maskb)
    
    end subroutine assign_vdw