program rmsd_main
use, intrinsic :: iso_fortran_env, only : output_unit, error_unit
use mctc_env, only : wp, error_type, fatal_error
use mctc_io, only : structure_type, read_structure
use targ
use rmsd
use rmsd_toml
use rmsd_config
use rmsd_filter
implicit none
character(len=*), parameter :: prog_name = "mctc-rmsd"
interface read_structure
procedure :: read_structure_from_arg
end interface read_structure
integer :: iarg, stat
type(arg_type), allocatable :: args(:)
type(error_type), allocatable :: error
type(structure_type) :: ref, mol
type(rmsd_config_type) :: config
type(rmsd_filter_type), pointer :: filter
type(toml_table), allocatable :: table, opts
type(toml_table), pointer :: child
type(toml_serializer) :: ser
character(len=:), allocatable :: rcfile, filter_name
real(wp) :: rmsd_val
type(targ_type) :: argparse
logical :: show_version, show_help, show_rc
logical, allocatable :: mask(:)
argparse = new_argument_parser()
call argparse%add_option("help")
call argparse%add_option("version")
call argparse%add_option("rc")
call argparse%add_option("filter", require=1)
call get_arguments(args)
call get_options(argparse, args, opts)
call get_value(opts, "filter", filter_name)
call get_value(opts, "version", show_version, .false.)
call get_value(opts, "rc", show_rc, .false.)
call get_value(opts, "help", show_help, size(args) <= 1 .and. .not.show_rc)
call opts%destroy
if (show_version) then
call version(output_unit)
stop
end if
if (show_help) then
call help(output_unit)
if (size(args) <= 1) error stop
stop
end if
call get_config_file(rcfile)
if (allocated(rcfile)) then
call read_config_file(table, rcfile, error)
if (allocated(error)) then
write(error_unit, '(a)') error%message
error stop
end if
else
table = toml_table()
end if
call get_value(table, "rmsd", child, requested=.true.)
if (associated(child)) then
call new_rmsd_config(config, child, error)
else
call fatal_error(error, "Type mismatch, rmsd must be a table")
end if
if (allocated(error)) then
write(error_unit, '(a)') error%message
error stop
end if
if (show_rc) then
if (allocated(rcfile)) then
write(output_unit, '(*(a:, 1x))') "#", "configuration file:", rcfile
else
write(output_unit, '(*(a:, 1x))') "#", "internal defaults"
end if
ser = toml_serializer(output_unit)
call table%accept(ser)
stop
end if
nullify(filter)
if (allocated(filter_name)) then
call config%get_filter(filter_name, filter)
if (.not.associated(filter)) then
write(error_unit, '(a)') &
& "No filter with name '"//filter_name//"' present"
error stop
end if
end if
call read_structure(ref, args(1), error)
if (allocated(error)) then
write(error_unit, '(a)') error%message
error stop
end if
if (associated(filter)) then
call filter%get_mask(ref, mask)
end if
stat = 0
do iarg = 2, size(args)
call read_structure(mol, args(iarg), error)
if (allocated(error)) then
write(error_unit, '(a)') error%message
stat = stat + 1
cycle
end if
if (config%strict) then
if (ref%nat /= mol%nat) then
write(error_unit, '(a)') &
"Number of atoms for '"//args(iarg)%arg//"' mismatch with reference"
stat = stat + 1
cycle
end if
if (any(ref%num(ref%id) /= mol%num(mol%id))) then
write(error_unit, '(a)') &
"Atomic number of '"//args(iarg)%arg//"' mismatch with reference"
stat = stat + 1
cycle
end if
end if
if (allocated(error)) then
stat = stat + 1
cycle
end if
call get_rmsd(ref, mol, rmsd_val, mask=mask)
write(output_unit, '(a, t40, es20.10, 1x, a)') &
& args(iarg)%arg, rmsd_val*config%conv, config%length_unit
end do
if (stat > 0) then
error stop
end if
contains
subroutine help(unit)
integer, intent(in) :: unit
write(unit, '(a, *(1x, a))') &
"Usage: "//prog_name//" [options] <file> <file>..."
write(unit, '(a)') &
"", &
"Requires at least two input files. Supported formats are:", &
"- xyz, mol, sdf, coord (0D), gen (C), pdb", &
"", &
"Configuration data is read from [rmsd] table in mctc.toml.", &
"Just place the configuration file mctc.toml (or .mctc.toml) in your home directory.", &
"Example:", &
"", &
" [rmsd]", &
" unit = ""AA""", &
" [rmsd.filter]", &
" heavy.exclude = [ ""H"", ""h"" ]", &
"", &
"Options", &
"-------", &
""
write(unit, '(3x, a, t25, a)') &
"--filter <name>", "Use <name> filter from configuration data to apply mask", &
"--rc", "check configuration data and print it to standard out", &
"--version", "Print program version and exit", &
"--help", "Show this help message"
write(unit, '(a)') &
"", &
"Filter", &
"------", &
"", &
"Filters can be defined in the [rmsd.filter] section, they take a list of", &
"atomic numbers and/or element symbols to define the allow-/deny-list.", &
"For example, to only check all carbon, nitrogen and oxygen atoms create", &
"a filter named organic with:", &
"", &
" organic.include = [6, 7, 8]", &
"", &
"Similarly, to create a filter for all heavy elements, effectively just", &
"excluding hydrogen with standard symbols, use:", &
"", &
" heavy.exclude = [""H"", ""h""]", &
"", &
"Note that this approach will still consider deuterium labeled as D,", &
"which would be excluded as well when using the atomic number instead.", &
"", &
"To create a PDB specific filter use the four character PDB identifier", &
"of the atoms and enable the PDB functionality.", &
"To match only the proteine backbone use", &
"", &
" c-alpha.include = ["" CA "", "" N "", "" C "", "" O ""]", &
" c-alpha.pdb = true", &
"", &
"Atomic numbers and element symbols can be included here as well.", &
""
end subroutine help
subroutine read_structure_from_arg(self, arg, error)
type(structure_type), intent(out) :: self
type(arg_type), intent(in) :: arg
type(error_type), allocatable, intent(out) :: error
call read_structure(self, arg%arg, error)
end subroutine read_structure_from_arg
subroutine get_variable(var, val)
!> Name of variable
character(len=*), intent(in) :: var
!> Value of variable
character(len=:), allocatable, intent(out) :: val
integer :: length, stat
call get_environment_variable(var, length=length, status=stat)
if (stat /= 0) then
return
endif
allocate(character(len=length) :: val, stat=stat)
if (stat /= 0) then
return
endif
if (length > 0) then
call get_environment_variable(var, val, status=stat)
if (stat /= 0) then
deallocate(val)
return
end if
end if
end subroutine get_variable
subroutine get_config_file(config)
!> Name of the configuration file
character(len=:), allocatable, intent(out) :: config
character(len=*), parameter :: rc = 'mctc.toml'
character(len=:), allocatable :: tmp, prefix
character :: sep
logical :: exist
if (is_windows()) then
sep = '\'
call get_variable('HOMEDRIVE', prefix)
call get_variable('HOMEDIR', tmp)
prefix = prefix // tmp
else
sep = '/'
call get_variable('HOME', prefix)
end if
if (allocated(prefix)) then
tmp = prefix // sep // rc
inquire(file=tmp, exist=exist)
if (exist) then
config = tmp
return
end if
tmp = prefix // sep // '.' // rc
inquire(file=tmp, exist=exist)
if (exist) then
config = tmp
return
end if
end if
inquire(file=rc, exist=exist)
if (exist) then
config = rc
return
end if
end subroutine get_config_file
!> Try to determine if we run on Windows and don't have POSIX compliance around
function is_windows() result(windows)
!> Operating system seems to be Windows
logical :: windows
character(len=:), allocatable :: tmp
windows = .false.
call get_variable('OS', tmp)
if (allocated(tmp)) then
windows = index(tmp, 'Windows_NT') > 0
end if
if (.not.windows) then
call get_variable('OSTYPE', tmp)
if (allocated(tmp)) then
windows = index(tmp, 'win') > 0 .or. index(tmp, 'msys') > 0
end if
end if
end function is_windows
subroutine version(unit)
integer, intent(in) :: unit
character(len=:), allocatable :: version_string
call get_rmsd_version(string=version_string)
write(unit, '(a, *(1x, a))') &
& prog_name, "version", version_string
end subroutine version
end program rmsd_main