| 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_urey(bds, prm_buf)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: urey_init
use mod_constants, only: angstrom2au, 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, iub, nub, &
cla, clb, clc, maxub, a, b, c, jc, jb
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classa(:), classb(:), classc(:), ubtmp(:)
real(rp), allocatable :: kub(:), l0ub(:)
logical :: 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
nub = 0
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:9) == 'ureybrad ') nub = nub + 1
end do
maxub = top%conn(2)%ri(top%mm_atoms+1)-1
! Maximum number of UB terms (each angle have an UB term)
call mallocate('assign_urey [classa]', nub, classa)
call mallocate('assign_urey [classb]', nub, classb)
call mallocate('assign_urey [classc]', nub, classc)
call mallocate('assign_urey [l0ub]', nub, l0ub)
call mallocate('assign_urey [kub]', nub, kub)
call mallocate('assign_urey [ubtmp]', maxub, ubtmp)
! Restart the reading from the beginning to actually save the parameters
iub = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:11) == 'urey-cubic ') then
tokb = 12
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong UREY-CUBIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%urey_cubic
! This parameter is 1/Angstrom
bds%urey_cubic = bds%urey_cubic / angstrom2au
else if(line(:13) == 'urey-quartic ') then
tokb = 14
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong UREY-QUARTIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%urey_quartic
bds%urey_quartic = bds%urey_quartic / (angstrom2au**2)
else if(line(:9) == 'ureybrad ') then
tokb = 10
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong UREYBRAD card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(iub)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong UREYBRAD card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(iub)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong UREYBRAD card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classc(iub)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong UREYBRAD card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kub(iub)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong UREYBRAD card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) l0ub(iub)
iub = iub + 1
end if
i = i+1
end do
ubtmp = -1
do a=1, top%mm_atoms
cla = top%atclass(a)
do jb=top%conn(2)%ri(a), top%conn(2)%ri(a+1)-1
b = top%conn(2)%ci(jb)
done = .false.
if(a > b) cycle
clb = top%atclass(b)
do jc=top%conn(1)%ri(a), top%conn(1)%ri(a+1)-1
c = top%conn(1)%ci(jc)
if(all(top%conn(1)%ci(top%conn(1)%ri(b): &
top%conn(1)%ri(b+1)-1) /= c)) cycle
! There is an angle in the form A-C-B
clc = top%atclass(c)
do j=1, nub
if((cla == classa(j) &
.and. clb == classc(j) &
.and. clc == classb(j)) .or. &
(clb == classa(j) &
.and. cla == classc(j) &
.and. clc == classb(j))) then
ubtmp(jb) = j
! Temporary assignament in a sparse matrix logic
done = .true.
exit
end if
end do
if(done) exit
! If we have already found a parameter for A-B pair, stop
! the research
end do
end do
end do
call urey_init(bds, count(ubtmp > 0))
iub = 1
do a=1, top%mm_atoms
do jb=top%conn(2)%ri(a), top%conn(2)%ri(a+1)-1
if(ubtmp(jb) > 0) then
bds%ureyat(1,iub) = a
bds%ureyat(2,iub) = top%conn(2)%ci(jb)
bds%kurey(iub) = kub(ubtmp(jb)) * kcalmol2au / (angstrom2au**2)
bds%l0urey(iub) = l0ub(ubtmp(jb)) * angstrom2au
iub = iub + 1
end if
end do
end do
call mfree('assign_urey [classa]', classa)
call mfree('assign_urey [classb]', classb)
call mfree('assign_urey [classc]', classc)
call mfree('assign_urey [l0ub]', l0ub)
call mfree('assign_urey [kub]', kub)
call mfree('assign_urey [ubtmp]', ubtmp)
end subroutine assign_urey