| 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(:) |
subroutine assign_pitors(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: pitors_init
use mod_constants, only: 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) :: il, i, tokb, toke, ipitors, npitors, &
cla, clb, maxpi, a, b, c, jb, iprm
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classa(:), classb(:), tmpat(:,:)
real(rp), allocatable :: kpi(:), tmpk(:)
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
npitors = 1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:7) == 'pitors ') npitors = npitors + 1
end do
maxpi = top%mm_atoms
! TODO This is maybe excessive, all trivalent atomso should be enough
call mallocate('assign_pitors [classa]', npitors, classa)
call mallocate('assign_pitors [classb]', npitors, classb)
call mallocate('assign_pitors [kpi]', npitors, kpi)
call mallocate('assign_pitors [tmpat]', 6, maxpi, tmpat)
call mallocate('assign_pitors [tmpk]', maxpi, tmpk)
! Restart the reading from the beginning to actually save the parameters
ipitors = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:7) == 'pitors ') then
tokb = 8
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong PITORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(ipitors)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong PITORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(ipitors)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong PITORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kpi(ipitors)
ipitors = ipitors + 1
end if
i = i+1
end do
ipitors = 1
do a=1, top%mm_atoms
! Check if the center is trigonal
if(top%conn(1)%ri(a+1) - top%conn(1)%ri(a) /= 3) cycle
cla = top%atclass(a)
! Loop over the atoms connected to the trignonal center
do jb=top%conn(1)%ri(a), top%conn(1)%ri(a+1)-1
b = top%conn(1)%ci(jb)
! This avoid to compute the functions two times
if(a > b) cycle
!Check if the second center is trigonal
if(top%conn(1)%ri(b+1) - top%conn(1)%ri(b) /= 3) cycle
clb = top%atclass(b)
do iprm=1, npitors
if((cla == classa(iprm) .and. clb == classb(iprm)) .or. &
(clb == classa(iprm) .and. cla == classb(iprm))) then
! The parameter is the right one
! Save the atoms in the following way:
!
! 2 5 a => 1
! \ / b => 4
! 1 -- 4
! / \
! 3 6
tmpat(:,ipitors) = 0
tmpat(1,ipitors) = a
do i=top%conn(1)%ri(a), top%conn(1)%ri(a+1)-1
c = top%conn(1)%ci(i)
if(c /= b) then
if(tmpat(2,ipitors) == 0) then
tmpat(2,ipitors) = c
else
tmpat(3,ipitors) = c
end if
end if
end do
tmpat(4,ipitors) = b
do i=top%conn(1)%ri(b), top%conn(1)%ri(b+1)-1
c = top%conn(1)%ci(i)
if(c /= a) then
if(tmpat(5,ipitors) == 0) then
tmpat(5,ipitors) = c
else
tmpat(6,ipitors) = c
end if
end if
end do
tmpk(ipitors) = kpi(iprm)
ipitors = ipitors+1
exit
end if
end do
end do
end do
call pitors_init(bds, ipitors-1)
do i=1, ipitors-1
bds%kpitors(i) = tmpk(i) * kcalmol2au
bds%pitorsat(:,i) = tmpat(:,i)
end do
call mfree('assign_pitors [classa]', classa)
call mfree('assign_pitors [classb]', classb)
call mfree('assign_pitors [kpi]', kpi)
call mfree('assign_pitors [tmpat]', tmpat)
call mfree('assign_pitors [tmpk]', tmpk)
end subroutine assign_pitors