assign_bond Subroutine

public subroutine assign_bond(bds, prm_buf, exclude_list, nexc_in)

Uses

  • proc~~assign_bond~~UsesGraph proc~assign_bond assign_bond module~mod_bonded mod_bonded proc~assign_bond->module~mod_bonded module~mod_constants mod_constants proc~assign_bond->module~mod_constants module~mod_io mod_io proc~assign_bond->module~mod_io module~mod_memory mod_memory proc~assign_bond->module~mod_memory module~mod_bonded->module~mod_io module~mod_bonded->module~mod_memory module~mod_topology mod_topology module~mod_bonded->module~mod_topology 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_topology->module~mod_memory module~mod_adjacency_mat mod_adjacency_mat module~mod_topology->module~mod_adjacency_mat module~mod_adjacency_mat->module~mod_memory

If there are no bonds in the system just return, there is nothing to do.

Arguments

Type IntentOptional Attributes Name
type(ommp_bonded_type), intent(inout) :: bds

Bonded potential data structure
character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)
integer(kind=ip), intent(in), optional, dimension(:) :: exclude_list

List of atoms for which interactions should not be computed
integer(kind=ip), intent(in), optional :: nexc_in

Number of atom in excluded list needed to skip a parameter

Calls

proc~~assign_bond~~CallsGraph proc~assign_bond assign_bond proc~read_atom_cards read_atom_cards proc~assign_bond->proc~read_atom_cards interface~mallocate mallocate proc~assign_bond->interface~mallocate proc~bond_init bond_init proc~assign_bond->proc~bond_init proc~isreal isreal proc~assign_bond->proc~isreal proc~tokenize tokenize proc~assign_bond->proc~tokenize proc~fatal_error fatal_error proc~assign_bond->proc~fatal_error proc~isint isint proc~assign_bond->proc~isint proc~ommp_message ommp_message proc~assign_bond->proc~ommp_message interface~mfree mfree proc~assign_bond->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~r_alloc3 r_alloc3 interface~mallocate->proc~r_alloc3 proc~i_alloc2 i_alloc2 interface~mallocate->proc~i_alloc2 proc~i_alloc1 i_alloc1 interface~mallocate->proc~i_alloc1 proc~r_alloc1 r_alloc1 interface~mallocate->proc~r_alloc1 proc~r_alloc2 r_alloc2 interface~mallocate->proc~r_alloc2 proc~i_alloc3 i_alloc3 interface~mallocate->proc~i_alloc3 proc~l_alloc1 l_alloc1 interface~mallocate->proc~l_alloc1 proc~l_alloc2 l_alloc2 interface~mallocate->proc~l_alloc2 proc~bond_init->interface~mallocate proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output proc~i_free1 i_free1 interface~mfree->proc~i_free1 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_free1 r_free1 interface~mfree->proc~r_free1 proc~r_free2 r_free2 interface~mfree->proc~r_free2 proc~l_free2 l_free2 interface~mfree->proc~l_free2 proc~i_free3 i_free3 interface~mfree->proc~i_free3 proc~chk_free chk_free proc~i_free1->proc~chk_free proc~l_free1->proc~chk_free proc~i_free2->proc~chk_free proc~r_free3->proc~chk_free proc~r_free1->proc~chk_free proc~chk_alloc chk_alloc proc~r_alloc3->proc~chk_alloc proc~memory_init memory_init proc~r_alloc3->proc~memory_init proc~i_alloc2->proc~chk_alloc proc~i_alloc2->proc~memory_init proc~i_alloc1->proc~chk_alloc proc~i_alloc1->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~i_alloc3->proc~chk_alloc proc~i_alloc3->proc~memory_init proc~l_alloc1->proc~chk_alloc proc~l_alloc1->proc~memory_init proc~l_alloc2->proc~chk_alloc proc~l_alloc2->proc~memory_init proc~close_output->proc~ommp_message proc~r_free2->proc~chk_free proc~l_free2->proc~chk_free proc~i_free3->proc~chk_free proc~chk_free->proc~fatal_error proc~chk_alloc->proc~fatal_error

Called by

proc~~assign_bond~~CalledByGraph proc~assign_bond assign_bond proc~init_bonded_for_link_atom init_bonded_for_link_atom proc~init_bonded_for_link_atom->proc~assign_bond proc~mmpol_init_from_xyz mmpol_init_from_xyz proc~mmpol_init_from_xyz->proc~assign_bond proc~ommp_system_from_qm_helper ommp_system_from_qm_helper proc~ommp_system_from_qm_helper->proc~assign_bond proc~ommp_create_link_atom ommp_create_link_atom proc~ommp_create_link_atom->proc~init_bonded_for_link_atom 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~ommp_init_xyz ommp_init_xyz proc~ommp_init_xyz->proc~mmpol_init_from_xyz proc~c_ommp_create_link_atom C_ommp_create_link_atom proc~c_ommp_create_link_atom->proc~ommp_create_link_atom proc~c_ommp_init_xyz C_ommp_init_xyz proc~c_ommp_init_xyz->proc~ommp_init_xyz

Contents

Source Code

assign_bond

Source Code

    subroutine assign_bond(bds, prm_buf, exclude_list, nexc_in)
        use mod_memory, only: mallocate, mfree
        use mod_io, only: fatal_error
        use mod_bonded, only: bond_init, ommp_bonded_type
        use mod_constants, only: angstrom2au, kcalmol2au
        
        implicit none

        type(ommp_bonded_type), intent(inout) :: bds
        !! Bonded potential data structure
        character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)
        integer(ip), dimension(:), intent(in), optional :: exclude_list
        !! List of atoms for which interactions should not be computed
        integer(ip), intent(in), optional :: nexc_in
        !! Number of atom in excluded list needed to skip a parameter

        integer(ip), parameter ::  nexc_default = 2
        integer(ip) :: il, i, j, l, jat, tokb, toke, ibnd, nbnd, &
                       cla, clb, nexc, iexc
        character(len=OMMP_STR_CHAR_MAX) :: line, errstring
        integer(ip), allocatable :: classa(:), classb(:)
        real(rp), allocatable :: kbnd(:), l0bnd(:)
        logical :: done
        type(ommp_topology_type), pointer :: top

        top => bds%top
        
        nexc = nexc_default
        if(present(exclude_list)) then
            if(present(nexc_in)) then
                if(nexc_in > 0 .and. nexc_in < 4) then
                    nexc = nexc_in
                end if
            end if
        end if

        if(.not. top%atclass_initialized .or. .not. top%atz_initialized) then
            call read_atom_cards(top, prm_buf)
        end if
        
        ! We assume that all pair of bonded atoms have a bonded 
        ! parameter
        call bond_init(bds, (top%conn(1)%ri(top%mm_atoms+1)-1) / 2)
        !! If there are no bonds in the system just return, there
        !! is nothing to do.
        if(.not. bds%use_bond) return

        bds%kbond = 0
        bds%l0bond = 0

        l=1
        do i=1, top%mm_atoms
            do j=top%conn(1)%ri(i), top%conn(1)%ri(i+1)-1
                jat = top%conn(1)%ci(j)
                if(i < jat) then
                    bds%bondat(1,l) = i
                    bds%bondat(2,l) = jat
                    l = l+1
                end if
            end do
        end do

        ! Read all the lines of file just to count how large vector should be 
        ! allocated 
        nbnd = 0
        do il=1, size(prm_buf) 
            line = prm_buf(il)
            if(line(:5) == 'bond ') nbnd = nbnd + 1
        end do

        call mallocate('assign_bond [classa]', nbnd, classa)
        call mallocate('assign_bond [classb]', nbnd, classb)
        call mallocate('assign_bond [l0bnd]', nbnd, l0bnd)
        call mallocate('assign_bond [kbnd]', nbnd, kbnd)

        ! Restart the reading from the beginning to actually save the parameters
        ibnd = 1
        i=1
        do il=1, size(prm_buf) 
            line = prm_buf(il)
           
            if(line(:11) == 'bond-cubic ') then
                tokb = 12
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong BOND-CUBIC card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) bds%bond_cubic
                ! This parameter is 1/Angstrom
                bds%bond_cubic = bds%bond_cubic / angstrom2au
            
            else if(line(:13) == 'bond-quartic ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong BOND-QUARTIC card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) bds%bond_quartic
                bds%bond_quartic = bds%bond_quartic / (angstrom2au**2)
            
            else if(line(:5) == 'bond ') then
                tokb = 6
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong BOND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) classa(ibnd)

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

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

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong BOND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) l0bnd(ibnd)
                
                ibnd = ibnd + 1
            end if
            i = i+1
        end do
        
        do i=1, size(bds%bondat,2)
            ! Atom class for current pair
            cla = top%atclass(bds%bondat(1,i))
            clb = top%atclass(bds%bondat(2,i))
            
            done = .false.
            do j=1, nbnd
                if((classa(j)==cla .and. classb(j)==clb) .or. &
                   (classa(j)==clb .and. classb(j)==cla)) then
                    done = .true.
                    bds%kbond(i) = kbnd(j) * kcalmol2au / (angstrom2au**2)
                    bds%l0bond(i) = l0bnd(j) * angstrom2au
                    exit
                end if
            end do

            if(present(exclude_list) .and. .not. done) then
                iexc = 0
                if(any(exclude_list == bds%bondat(1,i))) iexc = iexc + 1
                if(any(exclude_list == bds%bondat(2,i))) iexc = iexc + 1
                if(iexc >= nexc) then
                    done = .true.
                    write(errstring, '(A,I0,A,I0,A,I0,A)') &
                        "Bond ", bds%bondat(1,i), '-', bds%bondat(2,i), &
                        " not found and ignored because ", iexc, " atoms are &
                        &in excluded list."
                    call ommp_message(errstring, OMMP_VERBOSE_HIGH)
                end if
            end if
            
            if(.not. done) then
                call fatal_error("Bond parameter not found!")
            end if
        end do
        
        call mfree('assign_bond [classa]', classa)
        call mfree('assign_bond [classb]', classb)
        call mfree('assign_bond [l0bnd]', l0bnd)
        call mfree('assign_bond [kbnd]', kbnd)
    
    end subroutine assign_bond