If there are no bonds in the system just return, there is nothing to do.
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(:) | |||
integer(kind=ip), | intent(in), | optional, | dimension(:) | :: | exclude_list |
List of atoms for which interactions should not be computed |
integer(kind=ip), | intent(in), | optional | :: | nexc_in |
Number of atom in excluded list needed to skip a parameter |
subroutine assign_bond(bds, prm_buf, exclude_list, nexc_in)
use mod_memory, only: mallocate, mfree
use mod_io, only: fatal_error
use mod_bonded, only: bond_init, ommp_bonded_type
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), dimension(:), intent(in), optional :: exclude_list
!! List of atoms for which interactions should not be computed
integer(ip), intent(in), optional :: nexc_in
!! Number of atom in excluded list needed to skip a parameter
integer(ip), parameter :: nexc_default = 2
integer(ip) :: il, i, j, l, jat, tokb, toke, ibnd, nbnd, &
cla, clb, nexc, iexc
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classa(:), classb(:)
real(rp), allocatable :: kbnd(:), l0bnd(:)
logical :: done
type(ommp_topology_type), pointer :: top
top => bds%top
nexc = nexc_default
if(present(exclude_list)) then
if(present(nexc_in)) then
if(nexc_in > 0 .and. nexc_in < 4) then
nexc = nexc_in
end if
end if
end if
if(.not. top%atclass_initialized .or. .not. top%atz_initialized) then
call read_atom_cards(top, prm_buf)
end if
! We assume that all pair of bonded atoms have a bonded
! parameter
call bond_init(bds, (top%conn(1)%ri(top%mm_atoms+1)-1) / 2)
!! If there are no bonds in the system just return, there
!! is nothing to do.
if(.not. bds%use_bond) return
bds%kbond = 0
bds%l0bond = 0
l=1
do i=1, top%mm_atoms
do j=top%conn(1)%ri(i), top%conn(1)%ri(i+1)-1
jat = top%conn(1)%ci(j)
if(i < jat) then
bds%bondat(1,l) = i
bds%bondat(2,l) = jat
l = l+1
end if
end do
end do
! Read all the lines of file just to count how large vector should be
! allocated
nbnd = 0
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:5) == 'bond ') nbnd = nbnd + 1
end do
call mallocate('assign_bond [classa]', nbnd, classa)
call mallocate('assign_bond [classb]', nbnd, classb)
call mallocate('assign_bond [l0bnd]', nbnd, l0bnd)
call mallocate('assign_bond [kbnd]', nbnd, kbnd)
! Restart the reading from the beginning to actually save the parameters
ibnd = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:11) == 'bond-cubic ') then
tokb = 12
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong BOND-CUBIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%bond_cubic
! This parameter is 1/Angstrom
bds%bond_cubic = bds%bond_cubic / angstrom2au
else if(line(:13) == 'bond-quartic ') then
tokb = 14
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong BOND-QUARTIC card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%bond_quartic
bds%bond_quartic = bds%bond_quartic / (angstrom2au**2)
else if(line(:5) == 'bond ') then
tokb = 6
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong BOND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(ibnd)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong BOND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(ibnd)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong BOND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kbnd(ibnd)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong BOND card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) l0bnd(ibnd)
ibnd = ibnd + 1
end if
i = i+1
end do
do i=1, size(bds%bondat,2)
! Atom class for current pair
cla = top%atclass(bds%bondat(1,i))
clb = top%atclass(bds%bondat(2,i))
done = .false.
do j=1, nbnd
if((classa(j)==cla .and. classb(j)==clb) .or. &
(classa(j)==clb .and. classb(j)==cla)) then
done = .true.
bds%kbond(i) = kbnd(j) * kcalmol2au / (angstrom2au**2)
bds%l0bond(i) = l0bnd(j) * angstrom2au
exit
end if
end do
if(present(exclude_list) .and. .not. done) then
iexc = 0
if(any(exclude_list == bds%bondat(1,i))) iexc = iexc + 1
if(any(exclude_list == bds%bondat(2,i))) iexc = iexc + 1
if(iexc >= nexc) then
done = .true.
write(errstring, '(A,I0,A,I0,A,I0,A)') &
"Bond ", bds%bondat(1,i), '-', bds%bondat(2,i), &
" not found and ignored because ", iexc, " atoms are &
&in excluded list."
call ommp_message(errstring, OMMP_VERBOSE_HIGH)
end if
end if
if(.not. done) then
call fatal_error("Bond parameter not found!")
end if
end do
call mfree('assign_bond [classa]', classa)
call mfree('assign_bond [classb]', classb)
call mfree('assign_bond [l0bnd]', l0bnd)
call mfree('assign_bond [kbnd]', kbnd)
end subroutine assign_bond