| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_electrostatics_type), | intent(inout), | target | :: | eel |
Electrostatics data structure to be initialized |
|
| character(len=OMMP_STR_CHAR_MAX), | intent(in) | :: | prm_buf(:) |
Char buffer containing the prm file loaded in RAM |
subroutine assign_pol(eel, prm_buf)
use mod_memory, only: mallocate, mfree, ip, rp
use mod_electrostatics, only: set_screening_parameters
use mod_constants, only: angstrom2au
implicit none
type(ommp_electrostatics_type), intent(inout), target :: eel
!! Electrostatics data structure 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, l, iat, tokb, toke, ipg
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: polat(:), pgspec(:,:)
real(rp), allocatable :: thf(:), isopol(:)
real(rp) :: usc(4), psc(4), pisc(4), dsc(4)
integer(ip) :: npolarize, ipolarize
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 polarization asignament.")
end if
! Read all the lines of file just to count how large vector should be
! allocated
npolarize = 0
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:9) == 'polarize ') npolarize = npolarize + 1
end do
call mallocate('read_prm [polat]', npolarize, polat)
call mallocate('read_prm [isopol]', npolarize, isopol)
call mallocate('read_prm [thf]', npolarize, thf)
call mallocate('read_prm [pgspec]', 8_ip, npolarize, pgspec)
pgspec = 0
ipolarize = 1
! 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) == 'polarization ') then
tokb = 14
toke = tokenize(line, tokb)
select case(line(tokb:toke))
case('mutual')
continue
case('direct')
call fatal_error("Polarization DIRECT is not supported")
case default
call fatal_error("Wrong POLARIZATION card")
end select
else if(line(:15) == 'polar-12-intra ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-12-INTRA card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) pisc(1)
else if(line(:15) == 'polar-13-intra ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-13-INTRA card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) pisc(2)
else if(line(:15) == 'polar-14-intra ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-14-INTRA card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) pisc(3)
else if(line(:15) == 'polar-15-intra ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-15-INTRA card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) pisc(4)
else if(line(:15) == 'polar-12-scale ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-12-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) psc(1)
else if(line(:15) == 'polar-13-scale ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-13-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) psc(2)
else if(line(:15) == 'polar-14-scale ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-14-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) psc(3)
else if(line(:15) == 'polar-15-scale ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLAR-15-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) psc(4)
else if(line(:16) == 'direct-11-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong DIRECT-11-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) dsc(1)
else if(line(:16) == 'direct-12-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong DIRECT-12-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) dsc(2)
else if(line(:16) == 'direct-13-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong DIRECT-13-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) dsc(3)
else if(line(:16) == 'direct-14-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong DIRECT-14-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) dsc(4)
else if(line(:16) == 'mutual-11-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong MUTUAL-11-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) usc(1)
else if(line(:16) == 'mutual-12-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong MUTUAL-12-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) usc(2)
else if(line(:16) == 'mutual-13-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong MUTUAL-13-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) usc(3)
else if(line(:16) == 'mutual-14-scale ') then
tokb = 17
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong MUTUAL-14-SCALE card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) usc(4)
else if(line(:9) == 'polarize ') then
tokb = 10 ! len of keyword + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong POLARIZE card"
call fatal_error(errstring)
endif
read(line(tokb:toke), *) polat(ipolarize)
if(polat(ipolarize) < 0) then
write(errstring, *) "Wrong POLARIZE card (specific atom not supported)"
call fatal_error(errstring)
end if
! Isotropic polarizability
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLARIZE card"
call fatal_error(errstring)
endif
read(line(tokb:toke), *) isopol(ipolarize)
! Thole dumping factor
tokb = toke+1
toke = tokenize(line, tokb)
if(isint(line(tokb:toke))) then
! If this card is skipped then it is HIPPO
write(errstring, *) "HIPPO FF is not currently supported"
call fatal_error(errstring)
else if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong POLARIZE card"
call fatal_error(errstring)
endif
read(line(tokb:toke), *) thf(ipolarize)
! Polarization group information
tokb = toke+1
toke = tokenize(line, tokb)
if(isreal(line(tokb:toke)) .and. .not. isint(line(tokb:toke))) then
! If there is an additional real modifier then it is AMOEBA+
write(errstring, *) "AMOEBA+ FF is not currently supported"
call fatal_error(errstring)
end if
j = 1
do while(toke > 0)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong POLARIZE card"
call fatal_error(errstring)
endif
read(line(tokb:toke), *) pgspec(j,ipolarize)
tokb = toke+1
toke = tokenize(line, tokb)
j = j + 1
end do
ipolarize = ipolarize + 1
end if
i = i+1
end do
if(eel%amoeba) then
call set_screening_parameters(eel, eel%mscale, psc, dsc, usc, pisc)
eel%mmat_polgrp = 0
else
call set_screening_parameters(eel, eel%mscale, psc, dsc, usc)
end if
ipg = 0
! Now assign the parameters to the atoms
do i=1, size(top%attype)
! Polarization
do j=1, npolarize
if(polat(j) == top%attype(i)) then
eel%pol(i) = isopol(j) * angstrom2au**3
!TODO Thole factors.
! Assign a polgroup label to each atom
if(eel%mmat_polgrp(i) == 0) then
ipg = ipg+1
eel%mmat_polgrp(i) = ipg
end if
! loop over the atoms connected to ith atom
do k=top%conn(1)%ri(i), top%conn(1)%ri(i+1)-1
iat = top%conn(1)%ci(k)
if(any(top%attype(iat) == pgspec(:,j))) then
! The two atoms are in the same group
if(eel%mmat_polgrp(iat) == 0) then
eel%mmat_polgrp(iat) = eel%mmat_polgrp(i)
else if(eel%mmat_polgrp(iat) /= eel%mmat_polgrp(i)) then
! TODO This code have never been tested, as no
! suitable case have been found
do l=1, top%mm_atoms
if(eel%mmat_polgrp(l) == 0) then
continue
else if(eel%mmat_polgrp(l) == eel%mmat_polgrp(iat) &
.or. eel%mmat_polgrp(l) == eel%mmat_polgrp(i)) then
eel%mmat_polgrp(l) = min(eel%mmat_polgrp(iat), eel%mmat_polgrp(i))
else if(eel%mmat_polgrp(l) > max(eel%mmat_polgrp(iat),eel%mmat_polgrp(i))) then
eel%mmat_polgrp(l) = eel%mmat_polgrp(l) - 1
else
continue
end if
end do
end if
end if
end do
end if
end do
end do
call mfree('read_prm [polat]', polat)
call mfree('read_prm [isopol]', isopol)
call mfree('read_prm [thf]', thf)
call mfree('read_prm [pgspec]', pgspec)
end subroutine assign_pol