| 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_torsion(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: torsion_init
use mod_constants, only: kcalmol2au, deg2rad, eps_rp
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, ji, period
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classa(:), classb(:), classc(:), classd(:), &
t_n(:,:), tmpat(:,:), tmpprm(:), tmpbuf(:,:)
real(rp), allocatable :: t_amp(:,:), t_pha(:,:)
real(rp) :: amp, phase, torsion_unit = 1.0
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) == 'torsion ') nt = nt + 1
end do
maxt = top%conn(3)%ri(top%mm_atoms+1)-1
call mallocate('assign_torsion [classa]', nt, classa)
call mallocate('assign_torsion [classb]', nt, classb)
call mallocate('assign_torsion [classc]', nt, classc)
call mallocate('assign_torsion [classd]', nt, classd)
call mallocate('assign_torsion [t_amp]', 6_ip, nt, t_amp)
call mallocate('assign_torsion [t_pha]', 6_ip, nt, t_pha)
call mallocate('assign_torsion [t_n]', 6_ip, nt, t_n)
call mallocate('assign_torsion [tmpat]', 4_ip, maxt, tmpat)
call mallocate('assign_torsion [tmpprm]', maxt, tmpprm)
t_amp = 0.0
t_pha = 0.0
t_n = 1
! 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(:12) == 'torsionunit ') then
tokb = 13
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong TORSIONUNIT card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) torsion_unit
else if(line(:8) == 'torsion ') then
tokb = 9
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORSION 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 TORSION 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 TORSION 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 TORSION card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classd(it)
ji = 1
t_n(:,it) = -1
do j=1, 6
tokb = toke + 1
toke = tokenize(line, tokb)
if(toke < 0) exit
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong TORSION card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) amp
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong TORSION card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) phase
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong TORSION card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) period
if(abs(amp) > eps_rp) then
t_amp(ji, it) = amp
t_pha(ji, it) = phase
t_n(ji, it) = period
ji = ji + 1
end if
end do
if(j == 1) then
! No parameter found
write(errstring, *) "Wrong TORSION card"
call fatal_error(errstring)
end if
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
! Extrem check to avoid memory errors.
if(it > maxt) then
call mallocate('assign_torsion [tmpbuf]', 4, maxt, tmpbuf)
tmpbuf(:,:) = tmpat(:,:)
call mfree('assign_torsion [tmpat]', tmpat)
call mallocate('assign_torsion [tmpat]', 4, maxt+maxt+1, tmpat)
tmpat(:,1:maxt) = tmpbuf(:,:)
call mfree('assign_torsion [tmpbuf]', tmpbuf)
call mallocate('assign_torsion [tmpbuf]', 1, maxt, tmpbuf)
tmpbuf(1,:) = tmpprm(:)
call mfree('assign_torsion [tmpprm]', tmpprm)
call mallocate('assign_torsion [tmpprm]', maxt+maxt+1, tmpprm)
tmpprm(1:maxt) = tmpbuf(1,:)
call mfree('assign_torsion [tmpbuf]', tmpbuf)
maxt = 2*maxt + 1
end if
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 torsion_init(bds, it-1)
do i=1, it-1
bds%torsionat(:,i) = tmpat(:,i)
bds%torsamp(:,i) = t_amp(:,tmpprm(i)) * kcalmol2au * torsion_unit
bds%torsphase(:,i) = t_pha(:,tmpprm(i)) * deg2rad
bds%torsn(:,i) = t_n(:,tmpprm(i))
end do
call mfree('assign_torsion [classa]', classa)
call mfree('assign_torsion [classb]', classb)
call mfree('assign_torsion [classc]', classc)
call mfree('assign_torsion [classd]', classd)
call mfree('assign_torsion [t_amp]', t_amp)
call mfree('assign_torsion [t_pha]', t_pha)
call mfree('assign_torsion [t_n]', t_n)
call mfree('assign_torsion [tmpat]', tmpat)
call mfree('assign_torsion [tmpprm]', tmpprm)
end subroutine assign_torsion