| 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(:) |
Char buffer containing the prm file loaded in RAM |
||
| 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_angle(bds, prm_buf, exclude_list, nexc_in)
use mod_memory, only: mallocate, mfree
use mod_bonded, only: OMMP_ANG_SIMPLE, &
OMMP_ANG_H0, &
OMMP_ANG_H1, &
OMMP_ANG_H2, &
OMMP_ANG_INPLANE, &
OMMP_ANG_INPLANE_H0, &
OMMP_ANG_INPLANE_H1
use mod_bonded, only: angle_init
use mod_constants, only: kcalmol2au, rad2deg, deg2rad
implicit none
type(ommp_bonded_type), intent(inout) :: bds
!! Bonded potential data structure
character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)
!! Char buffer containing the prm file loaded in RAM
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 = 3
integer(ip) :: il, i, j, tokb, toke, iang, nang, &
cla, clb, clc, maxang, a, b, c, jc, jb, k, nhenv, &
iexc, nexc
character(len=OMMP_STR_CHAR_MAX) :: line, errstring
integer(ip), allocatable :: classa(:), classb(:), classc(:), angtype(:)
real(rp), allocatable :: kang(:), th0ang(:)
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
! Read all the lines of file just to count how large vector should be
! allocated
nang = 1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:6) == 'angle ') nang = nang + 3
! One angle keyourd could stand for 3 parameters for different H-env
if(line(:7) == 'anglep ') nang = nang + 2
! One angle keyourd could stand for 2 parameters for different H-env
end do
maxang = (top%conn(2)%ri(top%mm_atoms+1)-1) / 2
call mallocate('assign_angle [classa]', nang, classa)
call mallocate('assign_angle [classb]', nang, classb)
call mallocate('assign_angle [classc]', nang, classc)
call mallocate('assign_angle [eqang]', nang, th0ang)
call mallocate('assign_angle [kang]', nang, kang)
call mallocate('assign_angle [angtype]', nang, angtype)
call angle_init(bds, maxang)
! Restart the reading from the beginning to actually save the parameters
iang = 1
i=1
do il=1, size(prm_buf)
line = prm_buf(il)
if(line(:12) == 'angle-cubic ') then
tokb = 13
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE-CUBIC card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%angle_cubic
bds%angle_cubic = bds%angle_cubic * rad2deg
else if(line(:14) == 'angle-quartic ') then
tokb = 15
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE-QUARTIC card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%angle_quartic
bds%angle_quartic = bds%angle_quartic * rad2deg**2
else if(line(:13) == 'angle-pentic ') then
tokb = 13
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE-PENTIC card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%angle_pentic
bds%angle_pentic = bds%angle_pentic * rad2deg**3
else if(line(:13) == 'angle-sextic ') then
tokb = 14
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE-SEXTIC card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) bds%angle_sextic
bds%angle_sextic = bds%angle_sextic * rad2deg**4
else if(line(:6) == 'angle ') then
tokb = 7
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classc(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kang(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) th0ang(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(toke < 0) then
! Only one angle parameter is specified so it is good
! for all the H-envirnoment
angtype(iang) = OMMP_ANG_SIMPLE
iang = iang + 1
else
! Three equilibrim angles are specified for three different
! H-environment
angtype(iang) = OMMP_ANG_H0
iang = iang + 1
classa(iang) = classa(iang-1)
classb(iang) = classb(iang-1)
classc(iang) = classc(iang-1)
kang(iang) = kang(iang-1)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) th0ang(iang)
angtype(iang) = OMMP_ANG_H1
tokb = toke + 1
toke = tokenize(line, tokb)
if(toke > 0) then
iang = iang + 1
classa(iang) = classa(iang-1)
classb(iang) = classb(iang-1)
classc(iang) = classc(iang-1)
kang(iang) = kang(iang-1)
angtype(iang) = OMMP_ANG_H2
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLE card "
call fatal_error(errstring)
end if
read(line(tokb:toke), *) th0ang(iang)
end if
iang = iang + 1
end if
else if(line(:7) == 'anglep ') then
tokb = 8
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLEP card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classa(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLEP card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classb(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isint(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLEP card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) classc(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLEP card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) kang(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLEP card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) th0ang(iang)
tokb = toke + 1
toke = tokenize(line, tokb)
if(toke < 0) then
! Only one angle parameter is specified so it is good
! for all the H-envirnoment
angtype(iang) = OMMP_ANG_INPLANE
iang = iang + 1
else
! Three equilibrim angles are specified for three different
! H-environment
angtype(iang) = OMMP_ANG_INPLANE_H0
iang = iang + 1
classa(iang) = classa(iang-1)
classb(iang) = classb(iang-1)
classc(iang) = classc(iang-1)
kang(iang) = kang(iang-1)
if(.not. isreal(line(tokb:toke))) then
write(errstring, *) "Wrong ANGLEP card"
call fatal_error(errstring)
end if
read(line(tokb:toke), *) th0ang(iang)
angtype(iang) = OMMP_ANG_INPLANE_H1
iang = iang + 1
end if
end if
i = i+1
end do
nang = iang
iang = 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)
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)
done = .false.
do j=1, nang
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
if(angtype(j) == OMMP_ANG_SIMPLE .or. &
angtype(j) == OMMP_ANG_INPLANE) then
! For those types no check of the H
! environment is required
done = .true.
exit
else
! Check the H-environment
nhenv = 0
do k=top%conn(1)%ri(c), top%conn(1)%ri(c+1)-1
if(top%atz(top%conn(1)%ci(k)) == 1) &
nhenv = nhenv + 1
end do
if(top%atz(a) == 1) nhenv = nhenv-1
if(top%atz(b) == 1) nhenv = nhenv-1
if(nhenv == 0 .and. ( &
angtype(j) == OMMP_ANG_H0 .or. &
angtype(j) == OMMP_ANG_INPLANE_H0)) then
done = .true.
exit
else if(nhenv == 1 .and. ( &
angtype(j) == OMMP_ANG_H1 .or. &
angtype(j) == OMMP_ANG_INPLANE_H1)) then
done = .true.
exit
else if(nhenv == 2 .and. (&
angtype(j) == OMMP_ANG_H2)) then
done = .true.
exit
end if
end if
end if
end do
if(present(exclude_list) .and. .not. done) then
iexc = 0
if(any(exclude_list == a)) iexc = iexc + 1
if(any(exclude_list == b)) iexc = iexc + 1
if(any(exclude_list == c)) iexc = iexc + 1
if(iexc >= nexc) then
write(errstring, '(A,I0,A,I0,A,I0,A,I0,A)') &
"Angle ", a, '-', b, '-', c, " not found &
&and ignored because ", iexc, " atoms are &
&in excluded list."
call ommp_message(errstring, OMMP_VERBOSE_HIGH)
done = .true.
end if
end if
if(done) then
bds%angleat(1,iang) = a
bds%angleat(2,iang) = c
bds%angleat(3,iang) = b
bds%anglety(iang) = angtype(j)
bds%kangle(iang) = kang(j) * kcalmol2au
bds%eqangle(iang) = th0ang(j) * deg2rad
! Find the auxiliary atom for inplane angles
if(bds%anglety(iang) == OMMP_ANG_INPLANE .or. &
bds%anglety(iang) == OMMP_ANG_INPLANE_H0 .or. &
bds%anglety(iang) == OMMP_ANG_INPLANE_H1) then
! Find the auxiliary atom used to define the
! projection plane
if(bds%top%conn(1)%ri(bds%angleat(2,iang)+1) - &
bds%top%conn(1)%ri(bds%angleat(2,iang)) /= 3) &
then
call ommp_message("Angle IN-PLANE defined for a &
&non-trigonal center, this is only &
&acceptable if link-atoms should be &
&defined later.", OMMP_VERBOSE_LOW)
bds%kangle(iang) = 0.0
else
do j=bds%top%conn(1)%ri(bds%angleat(2,iang)), &
bds%top%conn(1)%ri(bds%angleat(2,iang)+1)-1
if(bds%top%conn(1)%ci(j) /= bds%angleat(1,iang) .and. &
bds%top%conn(1)%ci(j) /= bds%angleat(3,iang)) then
bds%angauxat(iang) = bds%top%conn(1)%ci(j)
endif
end do
end if
end if
iang = iang + 1
else
write(errstring, *) "No angle parameter found for &
&atoms ", a, b, c
call fatal_error(errstring)
end if
end do
end do
end do
call mfree('assign_angle [classa]', classa)
call mfree('assign_angle [classb]', classb)
call mfree('assign_angle [classc]', classc)
call mfree('assign_angle [eqang]', th0ang)
call mfree('assign_angle [kang]', kang)
call mfree('assign_angle [angtype]', angtype)
end subroutine assign_angle