assign_mpoles Subroutine

public subroutine assign_mpoles(eel, prm_buf)

Uses

  • proc~~assign_mpoles~~UsesGraph proc~assign_mpoles assign_mpoles module~mod_electrostatics mod_electrostatics proc~assign_mpoles->module~mod_electrostatics module~mod_constants mod_constants proc~assign_mpoles->module~mod_constants module~mod_memory mod_memory proc~assign_mpoles->module~mod_memory module~mod_electrostatics->module~mod_constants module~mod_electrostatics->module~mod_memory module~mod_profiling mod_profiling module~mod_electrostatics->module~mod_profiling module~mod_adjacency_mat mod_adjacency_mat module~mod_electrostatics->module~mod_adjacency_mat module~mod_topology mod_topology module~mod_electrostatics->module~mod_topology module~fmmlib_interface fmmlib_interface module~mod_electrostatics->module~fmmlib_interface module~mod_io mod_io module~mod_electrostatics->module~mod_io iso_c_binding iso_c_binding module~mod_constants->iso_c_binding module~mod_memory->module~mod_constants module~mod_memory->module~mod_io module~mod_memory->iso_c_binding module~mod_profiling->module~mod_constants module~mod_profiling->module~mod_memory module~mod_profiling->module~mod_io module~mod_adjacency_mat->module~mod_memory module~mod_topology->module~mod_memory module~mod_topology->module~mod_adjacency_mat module~fmmlib_interface->module~mod_constants module~mod_tree mod_tree module~fmmlib_interface->module~mod_tree module~mod_ribtree mod_ribtree module~fmmlib_interface->module~mod_ribtree module~mod_harmonics mod_harmonics module~fmmlib_interface->module~mod_harmonics module~mod_fmm_utils mod_fmm_utils module~fmmlib_interface->module~mod_fmm_utils module~mod_fmm mod_fmm module~fmmlib_interface->module~mod_fmm module~mod_octatree mod_octatree module~fmmlib_interface->module~mod_octatree module~mod_io->module~mod_constants module~mod_tree->module~mod_constants module~mod_tree->module~mod_adjacency_mat module~mod_tree->module~mod_fmm_utils module~mod_ribtree->module~mod_constants module~mod_ribtree->module~mod_profiling module~mod_ribtree->module~mod_tree module~mod_ribtree->module~mod_fmm_utils module~mod_harmonics->module~mod_constants module~mod_harmonics->module~mod_fmm_utils module~mod_fmm_utils->module~mod_constants module~mod_fmm->module~mod_constants module~mod_fmm->module~mod_tree module~mod_fmm->module~mod_harmonics module~mod_fmm->module~mod_fmm_utils module~mod_octatree->module~mod_constants module~mod_octatree->module~mod_profiling module~mod_octatree->module~mod_tree module~mod_octatree->module~mod_fmm_utils

Arguments

Type IntentOptional Attributes Name
type(ommp_electrostatics_type), intent(inout) :: eel

The electrostatic object to be initialized

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

Char buffer containing the prm file loaded in RAM


Calls

proc~~assign_mpoles~~CallsGraph proc~assign_mpoles assign_mpoles proc~fatal_error fatal_error proc~assign_mpoles->proc~fatal_error interface~mallocate mallocate proc~assign_mpoles->interface~mallocate proc~tokenize tokenize proc~assign_mpoles->proc~tokenize proc~isreal isreal proc~assign_mpoles->proc~isreal proc~isint isint proc~assign_mpoles->proc~isint proc~set_screening_parameters set_screening_parameters proc~assign_mpoles->proc~set_screening_parameters proc~ommp_message ommp_message proc~assign_mpoles->proc~ommp_message interface~mfree mfree proc~assign_mpoles->interface~mfree proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output 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~set_screening_parameters->proc~fatal_error 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~close_output->proc~ommp_message 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~r_free1->proc~chk_free 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_mpoles~~CalledByGraph proc~assign_mpoles assign_mpoles proc~init_eel_for_link_atom init_eel_for_link_atom proc~init_eel_for_link_atom->proc~assign_mpoles proc~mmpol_init_from_xyz mmpol_init_from_xyz proc~mmpol_init_from_xyz->proc~assign_mpoles proc~ommp_system_from_qm_helper ommp_system_from_qm_helper proc~ommp_system_from_qm_helper->proc~assign_mpoles proc~ommp_create_link_atom ommp_create_link_atom proc~ommp_create_link_atom->proc~init_eel_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


Source Code

    subroutine assign_mpoles(eel, prm_buf)
        use mod_memory, only: mallocate, mfree
        use mod_electrostatics, only: set_screening_parameters
        use mod_constants, only: AMOEBA_ROT_NONE, &
                                 AMOEBA_ROT_Z_THEN_X, &
                                 AMOEBA_ROT_BISECTOR, &
                                 AMOEBA_ROT_Z_ONLY, &
                                 AMOEBA_ROT_Z_BISECT, &
                                 AMOEBA_ROT_3_FOLD, &
                                 eps_rp, &
                                 OMMP_VERBOSE_DEBUG
        
        implicit none
        
        type(ommp_electrostatics_type), intent(inout) :: eel
        !! The electrostatic object to be initialized
        character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)
        !! Char buffer containing the prm file loaded in RAM

        integer(ip) :: il, i, j, k, iat, tokb, toke
        character(len=OMMP_STR_CHAR_MAX) :: line, errstring
        integer(ip), allocatable :: multat(:), multax(:,:), multframe(:)
        real(rp), allocatable :: cmult(:,:)
        real(rp) :: msc(4), csc(4) 
        real(rp) :: eel_au2kcalmol, eel_scale
        real(rp) :: default_eel_au2kcalmol = 332.063713
        ! Default conversion from A.U. to kcal/mol used in electrostatics of
        ! Tinker, only used to handle electric keyword
        integer(ip) :: nmult, nchg, imult, iax(3)
        logical :: ax_found(3), found13, only12, done
        type(ommp_topology_type), pointer :: top

        top => eel%top
        
        if(.not. top%attype_initialized) then
            call fatal_error("Atom type array in topology should be initialized&
                            & before performing multipoles asignament.")
        end if

        ! Read all the lines of file just to count how large vector should be 
        ! allocated 
        nmult = 0
        nchg = 0
        do il=1, size(prm_buf) 
            line = prm_buf(il)
            if(line(:11) == 'multipole ') nmult = nmult + 1
            if(line(:7) == 'charge ') nchg = nchg + 1
        end do
        
        ! MULTIPOLE
        call mallocate('read_prm [multat]', nmult+nchg, multat)
        call mallocate('read_prm [multframe]', nmult+nchg, multframe)
        call mallocate('read_prm [multax]', 3_ip, nmult+nchg, multax)
        call mallocate('read_prm [cmult]', 10_ip, nmult+nchg, cmult)
        multax = AMOEBA_ROT_NONE
        imult = 1
        
        ! Default values from Tinker manual
        msc = [0.0, 0.0, 1.0, 1.0]
        csc = [0.0, 0.0, 1.0, 1.0] 
        eel_scale = 1.0

        ! Restart the reading from the beginning to actually save the parameters
        do il=1, size(prm_buf) 
            line = prm_buf(il)
            
            if(line(:13) == 'chg-12-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong CHG-12-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) csc(1)
                if(csc(1) > 1.0) csc(1) = 1/csc(1)
            
            else if(line(:13) == 'chg-13-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong CHG-13-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) csc(2)
                if(csc(2) > 1.0) csc(2) = 1/csc(2)
            
            else if(line(:13) == 'chg-14-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong CHG-14-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) csc(3)
                if(csc(3) > 1.0) csc(3) = 1/csc(3)
            
            else if(line(:13) == 'chg-15-scale ') then
                tokb = 14
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong CHG-15-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) csc(4)
                if(csc(4) > 1.0) csc(4) = 1/csc(4)
            
            else if(line(:15) == 'mpole-12-scale ') then
                tokb = 16
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong MPOLE-12-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) msc(1)
                if(msc(1) > 1.0) msc(1) = 1/msc(1)
            
            else if(line(:15) == 'mpole-13-scale ') then
                tokb = 16
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong MPOLE-13-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) msc(2)
                if(msc(2) > 1.0) msc(2) = 1/msc(2)
            
            else if(line(:15) == 'mpole-14-scale ') then
                tokb = 16
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong MPOLE-14-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) msc(3)
                if(msc(3) > 1.0) msc(3) = 1/msc(3)
            
            else if(line(:15) == 'mpole-15-scale ') then
                tokb = 16
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong MPOLE-15-SCALE card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) msc(4)
                if(msc(4) > 1.0) msc(4) = 1/msc(4)
            
            else if(line(:9) == 'electric ') then
                tokb = 10
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong ELECTRIC card"
                    call fatal_error(errstring)
                end if
                ! This keyword is used to change the conversion from A.U. to
                ! kcal/mol for the electrostatic interaction.
                ! It is a bit creepy and unclear what it's correct to do here,
                ! I think that best thing is to scale the electrostatic itself
                ! by a factor (electric/default_electric) ** 0.5
                read(line(tokb:toke), *) eel_au2kcalmol
                eel_scale = (eel_au2kcalmol/default_eel_au2kcalmol) ** 0.5

            else if(line(:7) == 'charge ') then
                tokb = 8 ! len of keyword + 1
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong CHARGE card"
                    call fatal_error(errstring)
                endif
                read(line(tokb:toke), *) multat(imult)
                if(multat(imult) < 0) then
                    write(errstring, *) "Wrong CHARGE card (specific atom not supported)"
                    call fatal_error(errstring)
                end if
                ! Rotation axis information, charge does not contain any
                multax(:,imult) = 0
                multframe(imult) = AMOEBA_ROT_NONE

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

                read(line(tokb:toke), *) cmult(1, imult)
                cmult(2:10, imult) = 0.0 ! Fixed dipole and quadrupole are not present
                imult = imult + 1

            else if(line(:11) == 'multipole ') then
                tokb = 12 ! len of keyword + 1
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong MULTIPOLE card"
                    call fatal_error(errstring)
                endif
                read(line(tokb:toke), *) multat(imult)
                if(multat(imult) < 0) then
                    write(errstring, *) "Wrong MULTIPOLE card (specific atom not supported)"
                    call fatal_error(errstring)
                end if
                
                ! Rotation axis information
                tokb = toke+1
                toke = tokenize(line, tokb, 1)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong MULTIPOLE card at least two &
                                        &integers for axys specification expected."
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) multax(1,imult)
                
                tokb = toke+1
                toke = tokenize(line, tokb, 1)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong MULTIPOLE card at least two &
                                        &integers for axys specification expected."
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) multax(2,imult)

                ! For some centers also y axis is specified (integer) otherwise
                ! this is the zeroth-order multipole (charge)
                tokb = toke+1
                toke = tokenize(line, tokb)
                if(isint(line(tokb:toke))) then
                    ! Read the y axis
                    read(line(tokb:toke), *) multax(3,imult)
                    ! Get another token, this should be the charge.
                    tokb = toke+1 
                    toke = tokenize(line, tokb)
                end if

                ! The type of rotation is encoded in the sign/nullness of 
                ! of the axis specification
                if(multax(1,imult) == 0) then
                    multframe(imult) = AMOEBA_ROT_NONE
                else if(multax(2, imult) == 0) then
                    multframe(imult) = AMOEBA_ROT_Z_ONLY
                else if(multax(1,imult) < 0 .and. multax(2,imult) < 0 &
                        .and. multax(3,imult) < 0) then
                    multframe(imult) = AMOEBA_ROT_3_FOLD
                else if(multax(1,imult) < 0 .or. multax(2,imult) < 0) then
                    multframe(imult) = AMOEBA_ROT_BISECTOR
                else if(multax(2,imult) < 0 .and. multax(3,imult) < 0) then
                    multframe(imult) = AMOEBA_ROT_Z_BISECT
                else
                    multframe(imult) = AMOEBA_ROT_Z_THEN_X
                end if
                ! Remove the encoded information after saving it in a reasonable
                ! way
                multax(:,imult) = abs(multax(:,imult))

                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong MULTIPOLE card"
                    call fatal_error(errstring)
                end if

                read(line(tokb:toke), *) cmult(1, imult)

                read(prm_buf(il+1), *) cmult(2:4, imult)
                read(prm_buf(il+2), *) cmult(5, imult)
                read(prm_buf(il+3), *) cmult(6:7, imult)
                read(prm_buf(il+4), *) cmult(8:10, imult)
                !il = il+4
                
                imult = imult + 1
            end if
        end do
        
        if(nmult > 0 .and. nchg == 0) then
            call set_screening_parameters(eel, msc, eel%pscale, eel%dscale, &
                                          eel%uscale, eel%pscale_intra)
        else if(nmult == 0 .and. nchg > 0) then
            call set_screening_parameters(eel, csc, eel%pscale, eel%dscale, &
                                          eel%uscale)
        else if(nmult > 0 .and. nchg > 0) then
            write(errstring, *) "Unexpected FF with both multipoles and charges"
            call fatal_error(errstring)
        end if
        

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,only12,j,found13,ax_found,iax,iat,done) 
        do i=1, size(top%attype)
            if(eel%amoeba) then
                eel%mol_frame(i) = 0
                eel%ix(i) = 0
                eel%iy(i) = 0
                eel%iz(i) = 0
                eel%q0(:,i) = 0.0
            end if
            eel%q(:,i) = 0.0

            ! Flag to check assignament
            done = .false.

            ! Multipoles
            only12 = .false. ! Only search for params based on 12 connectivity
            do j=1, max(nmult, nchg)
                found13 = .false. ! Parameter found is based on 13 connectivity
                if(multat(j) == top%attype(i)) then
                    ! For each center different multipoles are defined for 
                    ! different environment. So first check if the environment
                    ! is the correct one
                    
                    ! Assignament with only 1,2-neighbours.
                    ax_found = .false.
                    iax = 0_ip

                    if(multframe(j) == AMOEBA_ROT_NONE) then
                        ! No axis needed
                        ax_found = .true.
                    else if(multframe(j) == AMOEBA_ROT_Z_ONLY) then
                        ! Assignament with only-z
                        ax_found(2:3) = .true.
                        do k=top%conn(1)%ri(i), top%conn(1)%ri(i+1)-1
                            iat = top%conn(1)%ci(k)
                            if(top%attype(iat) == multax(1,j) &
                               .and. .not. ax_found(1)) then
                                ax_found(1) = .true.
                                iax(1) = iat
                            end if
                        end do
                    else
                        ! 2 or 3 axis needed
                        if(multax(3,j) == 0) ax_found(3) = .true.
                        
                        ! Using only 1,2 connectivity
                        do k=top%conn(1)%ri(i), top%conn(1)%ri(i+1)-1
                            iat = top%conn(1)%ci(k)
                            if(top%attype(iat) == multax(1,j) &
                               .and. .not. ax_found(1)) then
                                ax_found(1) = .true.
                                iax(1) = iat
                            else if(top%attype(iat) == multax(2,j) &
                                    .and. .not. ax_found(2)) then
                                ax_found(2) = .true.
                                iax(2) = iat
                            else if(top%attype(iat) == multax(3,j) &
                                    .and. .not. ax_found(3)) then
                                ax_found(3) = .true.
                                iax(3) = iat
                            end if
                        end do

                        ! Using also 1,3 connectivity
                        if(ax_found(1) .and. .not. ax_found(2)) then
                            do k=top%conn(1)%ri(iax(1)), top%conn(1)%ri(iax(1)+1)-1
                                iat = top%conn(1)%ci(k)
                                if(iat == i .or. iat == iax(1)) cycle
                                if(top%attype(iat) == multax(2,j) &
                                   .and. .not. ax_found(2) &
                                   .and. iat /= iax(1)) then
                                    ax_found(2) = .true.
                                    iax(2) = iat
                                else if(top%attype(iat) == multax(3,j) &
                                        .and. .not. ax_found(3) & 
                                        .and. iat /= iax(1) &
                                        .and. iat /= iax(2)) then
                                    ax_found(3) = .true.
                                    iax(3) = iat
                                end if
                            end do
                            if(all(ax_found)) found13 = .true.
                        end if
                    end if

                    ! Everything is done, no further improvement is possible
                    if(all(ax_found) .and. .not. (only12 .and. found13)) then
                        if(eel%amoeba) then
                            eel%ix(i) = iax(2)
                            eel%iy(i) = iax(3)
                            eel%iz(i) = iax(1)
                            eel%mol_frame(i) = multframe(j)
                            eel%q(:,i) = cmult(:,j) 
                        else
                            eel%q(1,i) = cmult(1,j) 
                        end if
                        if(.not. done) then
                            done = .true.
                        else
                            write(errstring, "(A, I0)") &
                                "Reassigning multipoles for atom ", i
                            call ommp_message(errstring, OMMP_VERBOSE_DEBUG)
                        end if

                        write(errstring, "(A, I0, A, I0, A, I0, A, I0, A, I0, A)") &
                            "Atom ", i, " is assigned multipole set ", j, &
                            " axes [ ", iax(2), "-", iax(3), "-", iax(1), " ]"
                        call ommp_message(errstring, OMMP_VERBOSE_DEBUG)

                        
                        if(.not. found13) then
                            exit ! No further improvement is possible
                        else
                            only12 = .true.
                        end if
                    end if
                end if
            end do
            if(.not. done) then
                write(errstring, "(A, I0)") &
                    "No multipoles parameter found for atom ", i
                call ommp_message(errstring, OMMP_VERBOSE_LOW)
                call ommp_message("Previous error is only acceptable &
                                  &if link atom is used to fix those &
                                  &parameter later on", OMMP_VERBOSE_LOW)
            end if
        end do

        if(abs(eel_scale - 1.0) > eps_rp) then
            write(errstring, '(A, F10.6)') "Scaling charges by", eel_scale
            call ommp_message(errstring, OMMP_VERBOSE_LOW)
            eel%q = eel%q * eel_scale
        end if
        
        call mfree('read_prm [multat]', multat)
        call mfree('read_prm [multframe]', multframe)
        call mfree('read_prm [multax]', multax)
        call mfree('read_prm [cmult]', cmult)
    
    end subroutine assign_mpoles