| Type | Intent | Optional | 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 |
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, nmult+nchg, multax)
call mallocate('read_prm [cmult]', 10, 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 &
¶meter 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