| 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_strtor(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: strtor_init
use mod_constants, only: kcalmol2au, angstrom2au
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, j, tokb, toke, it, nt, &
cla, clb, clc, cld, maxt, a, b, c, d, jb, jc, jd, iprm
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classa(:), classb(:), classc(:), classd(:), &
tmpat(:,:), tmpprm(:)
real(rp), allocatable :: kat(:,:)
logical :: tor_done, bnd1_done, bnd2_done, bnd3_done
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
nt = 1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:8) == 'strtors ') nt = nt + 1
end do
maxt = top%conn(4)%ri(top%mm_atoms+1)-1
call mallocate('assign_strtor [classa]', nt, classa)
call mallocate('assign_strtor [classb]', nt, classb)
call mallocate('assign_strtor [classc]', nt, classc)
call mallocate('assign_strtor [classd]', nt, classd)
call mallocate('assign_strtor [kat]', 9_ip, nt, kat)
call mallocate('assign_strtor [tmpat]', 4_ip, maxt, tmpat)
call mallocate('assign_strtor [tmpprm]', maxt, tmpprm)
! Restart the reading from the beginning to actually save the parameters
it = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:8) == 'strtors ') then
tokb = 9
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong STRTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(it)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong STRTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(it)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong STRTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classc(it)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong STRTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classd(it)
do j=1, 9
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong STRTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kat(j,it)
end do
it = it + 1
end if
i = i+1
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
if(a > d) cycle
cld = top%atclass(d)
! There is a dihedral A-B-C-D
do iprm=1, nt
if((classa(iprm) == cla .and. &
classb(iprm) == clb .and. &
classc(iprm) == clc .and. &
classd(iprm) == cld) .or. &
(classa(iprm) == cld .and. &
classb(iprm) == clc .and. &
classc(iprm) == clb .and. &
classd(iprm) == cla)) then
! The parameter is ok
tmpat(:,it) = [a, b, c, d]
tmpprm(it) = iprm
it = it+1
exit
end if
end do
end do
end do
end do
end do
call strtor_init(bds, it-1)
do i=1, it-1
bds%strtorat(:,i) = tmpat(:,i)
if(classa(tmpprm(i)) == top%atclass(bds%strtorat(1,i))) then
bds%strtork(:,i) = kat(:,tmpprm(i))
else
bds%strtork(1:3,i) = kat(7:9,tmpprm(i))
bds%strtork(4:6,i) = kat(4:6,tmpprm(i))
bds%strtork(7:9,i) = kat(1:3,tmpprm(i))
end if
bds%strtork(:,i) = bds%strtork(:,i) * kcalmol2au / angstrom2au
tor_done = .false.
do j=1, size(bds%torsionat, 2)
if(all(bds%strtorat(:,i) == bds%torsionat(:,j))) then
tor_done = .true.
bds%strtor_t(i) = j
exit
end if
end do
bnd1_done = .false.
bnd2_done = .false.
bnd3_done = .false.
do j=1, size(bds%bondat, 2)
if(all(bds%strtorat(1:2,i) == bds%bondat(:,j)) .or. &
all(bds%strtorat(2:1:-1,i) == bds%bondat(:,j))) then
bnd1_done = .true.
bds%strtor_b(1,i) = j
else if(all(bds%strtorat(2:3,i) == bds%bondat(:,j)) .or. &
all(bds%strtorat(3:2:-1,i) == bds%bondat(:,j))) then
bnd2_done = .true.
bds%strtor_b(2,i) = j
else if(all(bds%strtorat(3:4,i) == bds%bondat(:,j)) .or. &
all(bds%strtorat(4:3:-1,i) == bds%bondat(:,j))) then
bnd3_done = .true.
bds%strtor_b(3,i) = j
end if
if(bnd1_done .and. bnd2_done .and. bnd3_done) exit
end do
if(.not. (tor_done .and. bnd1_done .and. bnd2_done .and. bnd3_done)) then
call fatal_error('Ill defined stretching-torsion coupling parameter')
end if
end do
call mfree('assign_strtor [classa]', classa)
call mfree('assign_strtor [classb]', classb)
call mfree('assign_strtor [classc]', classc)
call mfree('assign_strtor [classd]', classd)
call mfree('assign_strtor [kat]', kat)
call mfree('assign_strtor [tmpat]', tmpat)
call mfree('assign_strtor [tmpprm]', tmpprm)
end subroutine assign_strtor