| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_bonded_type), | intent(inout) | :: | bds |
Bonded potential data structure |
||
| character(len=OMMP_STR_CHAR_MAX), | intent(in) | :: | prm_buf(:) |
Char buffer containing the prm file loaded in RAM |
subroutine assign_tortors(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: tortor_newmap, tortor_init
use mod_constants, only: deg2rad, kcalmol2au
implicit none
type(ommp_bonded_type), intent(inout) :: bds
!! Bonded potential data 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, tokb, toke, iprm, jd, je, e, d, cle,it,cld,&
cla, clb, clc, a, b, c, jc, jb, itt, ndata, ntt, ibeg, iend, maxtt
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classx(:,:), map_dimension(:,:), tmpat(:,:), tmpprm(:), savedmap(:)
real(rp), allocatable :: data_map(:), ang_map(:,:)
type(ommp_topology_type), pointer :: top
top => bds%top
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
ntt = 0
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:8) == 'tortors ') ntt = ntt + 1
end do
maxtt = top%conn(4)%ri(top%mm_atoms+1)-1
call mallocate('assign_tortors [classx]', 5, ntt, classx)
call mallocate('assign_tortors [map_dimension]', 2, ntt, map_dimension)
call mallocate('assign_tortors [savedmap]', ntt, savedmap)
call mallocate('assign_tortors [tmpat]', 5, maxtt, tmpat)
call mallocate('assign_tortors [tmpprm]', maxtt, tmpprm)
! Restart the reading from the beginning to actually save the parameters
itt = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:8) == 'tortors ') then
tokb = 9
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classx(1,itt)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classx(2,itt)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classx(3,itt)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classx(4,itt)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classx(5,itt)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) map_dimension(1,itt)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) map_dimension(2,itt)
itt = itt + 1
end if
i = i+1
end do
! Allocate data space and finally read the map
ndata = dot_product(map_dimension(1,:), map_dimension(2,:))
call mallocate('assign_tortors [data_map]', ndata, data_map)
call mallocate('assign_tortors [ang_map]', 2, ndata, ang_map)
itt = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:8) == 'tortors ') then
ndata = map_dimension(1,itt)*map_dimension(2,itt)
do j=1, ndata
line = prm_buf(il+j)
tokb = tokenize(line)
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS data card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) ang_map(1, i)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS data card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) ang_map(2, i)
tokb = toke+1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong TORTORS data card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) data_map(i)
i = i + 1
end do
itt = itt + 1
end if
end do
it = 1
do a=1, top%mm_atoms
cla = top%atclass(a)
do jb=top%conn(1)%ri(a), top%conn(1)%ri(a+1)-1
b = top%conn(1)%ci(jb)
clb = top%atclass(b)
do jc=top%conn(1)%ri(b), top%conn(1)%ri(b+1)-1
c = top%conn(1)%ci(jc)
if(c == a) cycle
clc = top%atclass(c)
do jd=top%conn(1)%ri(c), top%conn(1)%ri(c+1)-1
d = top%conn(1)%ci(jd)
if(d == a .or. d == b) cycle
cld = top%atclass(d)
do je=top%conn(1)%ri(d), top%conn(1)%ri(d+1)-1
e = top%conn(1)%ci(je)
if(e == a .or. e == b .or. e == c) cycle
if(a > e) cycle
cle = top%atclass(e)
! There is a dihedral pair A-B-C-D-E
do iprm=1, ntt
if((classx(1,iprm) == cla .and. &
classx(2,iprm) == clb .and. &
classx(3,iprm) == clc .and. &
classx(4,iprm) == cld .and. &
classx(5,iprm) == cle) .or. &
(classx(1,iprm) == cle .and. &
classx(2,iprm) == cld .and. &
classx(3,iprm) == clc .and. &
classx(4,iprm) == clb .and. &
classx(5,iprm) == cla)) then
! The parameter is ok
tmpat(:,it) = [a, b, c, d, e]
tmpprm(it) = iprm
it = it+1
exit
end if
end do
end do
end do
end do
end do
end do
call tortor_init(bds, it-1)
savedmap = -1
iprm = 1
do i=1, it-1
if(savedmap(tmpprm(i)) < 0) then
! If needed, save the map in the module
ibeg = 1
do j=1, tmpprm(i)-1
ibeg = ibeg + map_dimension(1,j)*map_dimension(2,j)
end do
iend = ibeg + map_dimension(1,tmpprm(i))*map_dimension(2,tmpprm(i)) - 1
call tortor_newmap(bds, map_dimension(1,tmpprm(i)), &
map_dimension(2,tmpprm(i)), &
ang_map(1,ibeg:iend) * deg2rad, &
ang_map(2,ibeg:iend) * deg2rad, &
data_map(ibeg:iend) * kcalmol2au)
savedmap(tmpprm(i)) = iprm
iprm = iprm + 1
end if
bds%tortorat(:,i) = tmpat(:,i)
bds%tortorprm(i) = savedmap(tmpprm(i))
end do
call mfree('assign_tortors [classx]', classx)
call mfree('assign_tortors [map_dimension]', map_dimension)
call mfree('assign_tortors [savedmap]', savedmap)
call mfree('assign_tortors [data_map]', data_map)
call mfree('assign_tortors [ang_map]', ang_map)
call mfree('assign_tortors [tmpat]', tmpat)
call mfree('assign_tortors [tmpprm]', tmpprm)
end subroutine assign_tortors