| 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_opb(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: opb_init
use mod_constants, only: kcalmol2au, rad2deg
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, iopb, nopb, &
cla, clb, clc, cld, maxopb, a, b, c, d, jc, jb, iprm
character(len=OMMP_STR_CHAR_MAX) :: line, errstring, opb_type
integer(ip), allocatable :: classa(:), classb(:), classc(:), &
classd(:), tmpat(:,:)
real(rp), allocatable :: kopbend(:), 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
! Tinker manual default
opb_type = "w-d-c"
! Read all the lines of file just to count how large vector should be
! allocated
nopb = 0
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:7) == 'opbend ') nopb = nopb + 1
end do
if(nopb == 0) then
! If there are no OPB terms, just stop here.
return
end if
maxopb = top%mm_atoms*3
call mallocate('assign_opb [classa]', nopb, classa)
call mallocate('assign_opb [classb]', nopb, classb)
call mallocate('assign_opb [classc]', nopb, classc)
call mallocate('assign_opb [classd]', nopb, classd)
call mallocate('assign_opb [kopbend]', nopb, kopbend)
call mallocate('assign_opb [tmpat]', 4_ip, maxopb, tmpat)
call mallocate('assign_opb [tmpk]', maxopb, tmpk)
! Restart the reading from the beginning to actually save the parameters
iopb = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:11) == 'opbendtype ') then
tokb = 12
toke = tokenize(line, tokb)
read(line(tokb:toke), *) opb_type
else if(line(:13) == 'opbend-cubic ') then
tokb = 14
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND-CUBIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%opb_cubic
bds%opb_cubic = bds%opb_cubic * rad2deg
else if(line(:15) == 'opbend-quartic ') then
tokb = 16
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND-QUARTIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%opb_quartic
bds%opb_quartic = bds%opb_quartic * rad2deg**2
else if(line(:14) == 'opbend-pentic ') then
tokb = 15
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND-PENTIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%opb_pentic
bds%opb_pentic = bds%opb_pentic * rad2deg**3
else if(line(:14) == 'opbend-sextic ') then
tokb = 15
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND-SEXTIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%opb_sextic
bds%opb_sextic = bds%opb_sextic * rad2deg**4
else if(line(:7) == 'opbend ') then
tokb = 8
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(iopb)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(iopb)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classc(iopb)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classd(iopb)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong OPBEND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kopbend(iopb)
iopb = iopb + 1
end if
i = i+1
end do
iopb = 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)
clb = top%atclass(b)
c = -1
d = -1
clc = 0
cld = 0
do jc=top%conn(1)%ri(a), top%conn(1)%ri(a+1)-1
if(top%conn(1)%ci(jc) == b) cycle
if(c < 0) then
c = top%conn(1)%ci(jc)
clc = top%atclass(c)
else if(d < 0) then
d = top%conn(1)%ci(jc)
cld = top%atclass(d)
end if
end do
do iprm=1, nopb
if((classa(iprm) == clb) .and. &
(classb(iprm) == cla) .and. &
((classc(iprm) == clc .and. &
classd(iprm) == cld) .or. &
(classd(iprm) == clc .and. &
classc(iprm) == cld) .or. &
(classd(iprm) == 0 .and. &
(classc(iprm) == cld .or. classc(iprm) == clc)) .or. &
(classc(iprm) == 0 .and. &
(classd(iprm) == cld .or. classd(iprm) == clc)) .or. &
(classc(iprm) == 0 .or. classd(iprm) == 0))) then
! The parameter is ok
tmpat(1,iopb) = a
tmpat(2,iopb) = b
tmpat(3,iopb) = c
tmpat(4,iopb) = d
tmpk(iopb) = kopbend(iprm)
iopb = iopb + 1
exit
endif
end do
end do
end do
call opb_init(bds, iopb-1, trim(opb_type))
do i=1, iopb-1
bds%kopb(i) = tmpk(i) * kcalmol2au
bds%opbat(:,i) = tmpat(:,i)
end do
call mfree('assign_opb [classa]', classa)
call mfree('assign_opb [classb]', classb)
call mfree('assign_opb [classc]', classc)
call mfree('assign_opb [classd]', classd)
call mfree('assign_opb [kopbend]', kopbend)
call mfree('assign_opb [tmpat]', tmpat)
call mfree('assign_opb [tmpk]', tmpk)
end subroutine assign_opb