| 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_angtor(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: angtor_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, 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, ang1_done, ang2_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) == 'angtors ') nt = nt + 1
end do
maxt = top%conn(4)%ri(top%mm_atoms+1)-1
call mallocate('assign_angtor [classa]', nt, classa)
call mallocate('assign_angtor [classb]', nt, classb)
call mallocate('assign_angtor [classc]', nt, classc)
call mallocate('assign_angtor [classd]', nt, classd)
call mallocate('assign_angtor [kat]', 6, nt, kat)
call mallocate('assign_angtor [tmpat]', 4, maxt, tmpat)
call mallocate('assign_angtor [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) == 'angtors ') then
tokb = 9
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGTORS 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 ANGTORS 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 ANGOTORS 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 ANGTORS card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classd(it)
do j=1, 6
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGTORS 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 angtor_init(bds, it-1)
do i=1, it-1
bds%angtorat(:,i) = tmpat(:,i)
if(classa(tmpprm(i)) == top%atclass(bds%angtorat(1,i))) then
bds%angtork(:,i) = kat(:,tmpprm(i)) * kcalmol2au
else
bds%angtork(1:3,i) = kat(4:6,tmpprm(i)) * kcalmol2au
bds%angtork(4:6,i) = kat(1:3,tmpprm(i)) * kcalmol2au
end if
tor_done = .false.
do j=1, size(bds%torsionat, 2)
if(all(bds%angtorat(:,i) == bds%torsionat(:,j))) then
tor_done = .true.
bds%angtor_t(i) = j
exit
end if
end do
ang1_done = .false.
ang2_done = .false.
do j=1, size(bds%angleat, 2)
if(all(bds%angtorat(1:3,i) == bds%angleat(:,j)) .or. &
all(bds%angtorat(1:3,i) == bds%angleat(3:1:-1,j))) then
ang1_done = .true.
bds%angtor_a(1,i) = j
else if(all(bds%angtorat(2:4,i) == bds%angleat(:,j)) .or. &
all(bds%angtorat(2:4,i) == bds%angleat(3:1:-1,j))) then
ang2_done = .true.
bds%angtor_a(2,i) = j
end if
if(ang1_done .and. ang2_done) exit
end do
if(.not. (tor_done .and. ang1_done .and. ang2_done)) then
call fatal_error('Ill defined angle-torsion coupling parameter')
end if
end do
call mfree('assign_angtor [classa]', classa)
call mfree('assign_angtor [classb]', classb)
call mfree('assign_angtor [classc]', classc)
call mfree('assign_angtor [classd]', classd)
call mfree('assign_angtor [kat]', kat)
call mfree('assign_angtor [tmpat]', tmpat)
call mfree('assign_angtor [tmpprm]', tmpprm)
end subroutine assign_angtor