This function read a Tinker xyz and prm files and initialize all the quantities need to describe the environment within this library.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(ommp_system), | intent(inout), | target | :: | sys_obj |
The system object to be initialized |
|
character(len=*), | intent(in) | :: | xyz_file |
name of the input XYZ file |
||
character(len=*), | intent(in) | :: | prm_file |
name of the input PRM file |
subroutine mmpol_init_from_xyz(sys_obj, xyz_file, prm_file)
!! This function read a Tinker xyz and prm files
!! and initialize all the quantities need to describe the environment
!! within this library.
use mod_mmpol, only: mmpol_prepare, mmpol_init, mmpol_init_nonbonded, &
mmpol_init_bonded
use mod_topology, only: ommp_topology_type, check_conn_matrix
use mod_electrostatics, only: ommp_electrostatics_type
use mod_memory, only: ip, mfree, mallocate, memory_init
use mod_constants, only: angstrom2au, OMMP_VERBOSE_DEBUG
use mod_adjacency_mat, only: adj_mat_from_conn, yale_sparse, &
build_conn_upto_n
use mod_prm, only: check_keyword, assign_pol, assign_mpoles, &
assign_vdw, assign_bond, assign_angle, assign_urey, &
assign_strbnd, assign_opb, assign_pitors, &
assign_torsion, assign_tortors, assign_angtor, &
assign_strtor, assign_imptorsion, get_prm_ff_type
use mod_utils, only: starts_with_alpha, isreal, isint, &
str_to_lower, str_uncomment, &
tokenize_pure, atoi, atof
use mod_io, only: large_file_read
implicit none
type(ommp_system), intent(inout), target :: sys_obj
!! The system object to be initialized
character(len=*), intent(in) :: xyz_file
!! name of the input XYZ file
character(len=*), intent(in) :: prm_file
!! name of the input PRM file
integer(ip), parameter :: iof_xyzinp = 200, &
maxn12 = 8
integer(ip) :: my_mm_atoms, ist, i, j, atom_id, lc, tokx(2)
integer(ip), allocatable :: i12(:,:), attype(:)
character(len=OMMP_STR_CHAR_MAX) :: msg
character(len=OMMP_STR_CHAR_MAX), allocatable :: lines(:), prm_buf(:)
logical :: fex
type(yale_sparse) :: adj
type(ommp_topology_type), pointer :: top
type(ommp_electrostatics_type), pointer :: eel
call time_push()
call time_push()
! Check that files are present
inquire(file=prm_file(1:len(trim(prm_file))), exist=fex)
if(.not. fex) then
call fatal_error("Input file '"//&
&prm_file(1:len(trim(prm_file)))//"' does not exist.")
end if
inquire(file=xyz_file(1:len(trim(xyz_file))), exist=fex)
if(.not. fex) then
call fatal_error("Input file '"//&
&xyz_file(1:len(trim(xyz_file)))//"' does not exist.")
end if
write(msg, "(A)") "Reading XYZ file: "//xyz_file(1:len(trim(xyz_file)))
call ommp_message(msg, OMMP_VERBOSE_DEBUG)
! read tinker xyz file
call time_push()
call large_file_read(xyz_file, lines)
call time_pull('XYZ File reading')
call large_file_read(prm_file, prm_buf)
! Remove comments from prm file
!$omp parallel do
do i=1, size(prm_buf)
prm_buf(i) = str_to_lower(prm_buf(i))
prm_buf(i) = str_uncomment(prm_buf(i), '!')
end do
! First line contains as first word the number of atoms and
! then a comment that could be ignored.
read(lines(1), *) my_mm_atoms
! Initialize the mmpol module
! TODO I'm assuming that it is AMOEBA and fully polarizable
call mmpol_init(sys_obj, get_prm_ff_type(prm_buf), &
my_mm_atoms, my_mm_atoms)
! Those are just shortcut
top => sys_obj%top
eel => sys_obj%eel
! Temporary quantities that are only used during the initialization
call mallocate('mmpol_init_from_xyz [attype]', my_mm_atoms, attype)
call mallocate('mmpol_init_from_xyz [i12]', maxn12, my_mm_atoms, i12)
call time_push()
!$omp parallel do default(shared) private(i, lc, tokx) &
!$omp private(atom_id, j)
do i=1, my_mm_atoms
lc = i+1
! Initializations
attype(i) = 0_ip
i12(:,i) = 0_ip
eel%polar_mm(i) = i
eel%pol(i) = 0.0
! First token contains an atom ID. Only sequential numbering is
! currently supported.
!dir$ forceinline
tokx = tokenize_pure(lines(lc))
!dir$ forceinline
tokx = tokenize_pure(lines(lc), tokx(2))
atom_id = atoi(lines(lc)(tokx(1):tokx(2)))
if(atom_id /= i) then
call fatal_error('Non-sequential atom ids in xyz cannot be handled')
end if
! This token should contain an atom name, so it should start
! with a letter. If this is not true, an unexpected error should
! be raised.
!dir$ forceinline
tokx = tokenize_pure(lines(lc), tokx(2)+1)
if(isreal(lines(lc)(tokx(1):tokx(2)))) then
call fatal_error('Atom symbol missing (or completely numerical) or PBC string present in XYZ')
end if
! The remaining part contains cartesian coordinates, atom type
! and connectivity for the current atom.
do j=1, 3
! Read coordinates and convert to angstrom
!dir$ forceinline
tokx = tokenize_pure(lines(lc), tokx(2)+1)
top%cmm(j,i) = atof(lines(lc)(tokx(1):tokx(2))) * angstrom2au
end do
!dir$ forceinline
tokx = tokenize_pure(lines(lc), tokx(2)+1)
attype(i) = atoi(lines(lc)(tokx(1):tokx(2)))
do j=1, maxn12
!dir$ forceinline
tokx = tokenize_pure(lines(lc), tokx(2)+1)
if(tokx(2) < 0) exit
i12(j,i) = atoi(lines(lc)(tokx(1):tokx(2)))
end do
end do
deallocate(lines)
call time_pull('Interpreting XYZ')
! Writes the adjacency matrix in Yale sparse format in adj and then
! build the connectivity up to 4th order. This is needed here to be
! able to assign the parameters
call adj_mat_from_conn(i12, adj)
call build_conn_upto_n(adj, 4, top%conn, .false.)
call mfree('mmpol_init_from_xyz [i12]', i12)
if( .not. check_keyword(prm_buf)) then
call fatal_error("PRM file cannot be completely understood")
end if
call time_pull("XYZ reading and topology generation")
top%attype = attype
top%attype_initialized = .true.
call mfree('mmpol_init_from_xyz [attype]', attype)
call time_push()
call ommp_message("Assigning electrostatic parameters", OMMP_VERBOSE_DEBUG)
call assign_pol(eel, prm_buf)
call assign_mpoles(eel, prm_buf)
call ommp_message("Assigning non-bonded parameters", OMMP_VERBOSE_DEBUG)
call mmpol_init_nonbonded(sys_obj)
call assign_vdw(sys_obj%vdw, top, prm_buf)
call ommp_message("Assigning bonded parameters", OMMP_VERBOSE_DEBUG)
call mmpol_init_bonded(sys_obj)
call check_conn_matrix(sys_obj%top, 4)
call assign_bond(sys_obj%bds, prm_buf)
call assign_angle(sys_obj%bds, prm_buf)
call assign_urey(sys_obj%bds, prm_buf)
call assign_strbnd(sys_obj%bds, prm_buf)
call assign_opb(sys_obj%bds, prm_buf)
call assign_pitors(sys_obj%bds, prm_buf)
call assign_torsion(sys_obj%bds, prm_buf)
call assign_imptorsion(sys_obj%bds, prm_buf)
call assign_tortors(sys_obj%bds, prm_buf)
call assign_angtor(sys_obj%bds, prm_buf)
call assign_strtor(sys_obj%bds, prm_buf)
call time_pull('Total prm assignament')
deallocate(prm_buf)
call mmpol_prepare(sys_obj)
call time_pull('MMPol initialization from .xyz file')
end subroutine mmpol_init_from_xyz