comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Added CIF (Christallographic Information File) as possible format for inputting

Browse source 
coordinates. Added 1 regtest and a couple of CIF examples. Only the basics of
the CIF have been implemented (it may require improvements in the future).


svn-origin-rev: 8095
This commit is contained in:
Teodoro Laino 2008-12-31 08:38:47 +00:00
parent 9eae5c2555
commit 39d24ab3ce
10 changed files with 861 additions and 38 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/OBJECTDEFS
View file

@ -595,6 +595,7 @@ OBJECTS_GENERIC =\
timings_mp.o\
topology.o\
topology_amber.o\
topology_cif.o\
topology_coordinate_util.o\
topology_connectivity_util.o\
topology_constraint_util.o\

11
src/input_constants.F
View file

@ -268,11 +268,12 @@ MODULE input_constants
do_conn_g87=6,&
do_conn_g96=7,&
do_conn_amb7=8,&
do_coord_pdb=0,&
do_coord_xyz=1,&
do_coord_g96=2,&
do_coord_crd=3,&
do_coord_off=4
do_coord_off=0,&
do_coord_pdb=1,&
do_coord_xyz=2,&
do_coord_g96=3,&
do_coord_crd=4,&
do_coord_cif=5
INTEGER, PARAMETER, PUBLIC :: do_ff_undef=0,&
do_ff_quartic=1,&

38
src/input_cp2k_subsys.F
View file

@ -333,38 +333,43 @@ CONTAINS
"controls the printing of information in the topology settings", &
print_level=high_print_level,filename="__STD_OUT__",&
error=error)
CALL keyword_create(keyword,"cif_info",&
description="if the printkey is active prints information when parsing CIF files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"pdb_info",&
description="if the printkey is active prints information regarding PDB",&
description="if the printkey is active prints information when parsing PDB files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"xyz_info",&
description="if the printkey is active prints information regarding XYZ",&
description="if the printkey is active prints information when parsing XYZ files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"psf_info",&
description="if the printkey is active prints information regarding PSF",&
description="if the printkey is active prints information when parsing PSF files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"amber_info",&
description="if the printkey is active prints information regarding ABER topology file",&
description="if the printkey is active prints information when parsing ABER topology files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"g96_info",&
description="if the printkey is active prints information regarding G96",&
description="if the printkey is active prints information when parsing G96 files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"crd_info",&
description="if the printkey is active prints information regarding CRD",&
description="if the printkey is active prints information when parsing CRD files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword,"gtop_info",&
description="if the printkey is active prints information regarding GROMOS Topology",&
description="if the printkey is active prints information when parsing GROMOS topology files.",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(print_key,keyword,error=error)
CALL keyword_release(keyword,error=error)
@ -1301,14 +1306,17 @@ CONTAINS
CALL keyword_create(keyword, name="COORD_FILE_FORMAT",&
variants=s2a("COORDINATE"),&
description="Set up the way in which coordinates will be read.",&
usage="COORD_FILE_FORMAT (PDB|XYZ|G96|CRD|OFF)", &
enum_c_vals=s2a( "PDB","XYZ","G96","CRD","OFF"),&
enum_i_vals=(/do_coord_pdb, do_coord_xyz, do_coord_g96, do_coord_crd, do_coord_off/),&
enum_desc=s2a("Coordinates provided through a PDB file format",&
"Coordinates provided through an XYZ file format",&
"Coordinates provided through a GROMOS96 file format",&
"Coordinates provided through an AMBER file format",&
"Coordinates read in the &COORD section of the input file"),&
usage="COORD_FILE_FORMAT (OFF|PDB|XYZ|G96|CRD|CIF)", &
enum_c_vals=s2a( "OFF","PDB","XYZ","G96","CRD","CIF"),&
enum_i_vals=(/do_coord_off, do_coord_pdb, do_coord_xyz, do_coord_g96, do_coord_crd,&
do_coord_cif/),&
enum_desc=s2a(&
"Coordinates read in the &COORD section of the input file",&
"Coordinates provided through a PDB file format",&
"Coordinates provided through an XYZ file format",&
"Coordinates provided through a GROMOS96 file format",&
"Coordinates provided through an AMBER file format",&
"Coordinates provided through a CIF (Crystallographic Information File) file format"),&
default_i_val=do_coord_off,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)

31
src/topology.F
View file

@ -21,7 +21,8 @@ MODULE topology
USE input_constants, ONLY: &
do_conn_amb7, do_conn_g87, do_conn_g96, do_conn_generate, &
do_conn_mol_set, do_conn_off, do_conn_psf, do_conn_psf_u, &
do_coord_crd, do_coord_g96, do_coord_pdb, do_coord_xyz
do_coord_cif, do_coord_crd, do_coord_g96, do_coord_off, do_coord_pdb, &
do_coord_xyz
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
@ -42,6 +43,7 @@ MODULE topology
timestop
USE topology_amber, ONLY: read_connectivity_amber,&
read_coordinate_crd
USE topology_cif, ONLY: read_coordinate_cif
USE topology_connectivity_util, ONLY: topology_connectivity_pack
USE topology_constraint_util, ONLY: topology_constraint_pack
USE topology_coordinate_util, ONLY: topology_coordinate_pack
@ -510,30 +512,31 @@ CONTAINS
!-----------------------------------------------------------------------------
CALL timeset(routineN//"_READ_COORDINATE",handle2)
SELECT CASE (topology%coord_type)
CASE (do_coord_off)
! Do nothing.. we will parse later from the &COORD section..
CASE (do_coord_g96)
CALL read_coordinate_g96 (topology,para_env,subsys_section,error)
CALL read_atoms_input ( topology, overwrite=.TRUE.,&
subsys_section=subsys_section,save_mem=save_mem,error=error)
CASE (do_coord_crd)
CALL read_coordinate_crd (topology,para_env,subsys_section,error)
CALL read_atoms_input ( topology, overwrite=.TRUE.,&
subsys_section=subsys_section,save_mem=save_mem,error=error)
CASE (do_coord_pdb)
CALL read_coordinate_pdb (topology,para_env,subsys_section,error)
CALL read_atoms_input ( topology, overwrite=.TRUE.,&
subsys_section=subsys_section,save_mem=save_mem,error=error)
CASE (do_coord_xyz)
CALL read_coordinate_xyz (topology,para_env,subsys_section,error)
CALL read_atoms_input ( topology, overwrite=.TRUE.,&
subsys_section=subsys_section,save_mem=save_mem,error=error)
CASE (do_coord_cif)
CALL read_coordinate_cif (topology,para_env,subsys_section,error)
CASE DEFAULT
CALL read_atoms_input ( topology, subsys_section=subsys_section,&
save_mem=save_mem,error=error)
! We should never reach this point..
CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure)
END SELECT
! Parse &COORD section and in case overwrite..
CALL read_atoms_input ( topology, overwrite=(topology%coord_type/=do_coord_off),&
subsys_section=subsys_section,save_mem=save_mem,error=error)
CALL timestop(handle2)
IF (topology%natoms.LE.0) THEN
CALL stop_program ("coordinate_control","No atomic coordinates have been found")
END IF
! Check on atom numbers
CALL cp_assert(topology%natoms>0,cp_fatal_level,cp_assertion_failed,routineP,&
"No atomic coordinates have been found! "//&
CPSourceFileRef,&
error,failure)
CALL timestop(handle)
CALL cp_print_key_finished_output(iw,logger,subsys_section,&
"PRINT%TOPOLOGY_INFO",error=error)

419
src/topology_cif.F Normal file
View file

@ -0,0 +1,419 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2008 CP2K developers group !
!-----------------------------------------------------------------------------!
! *****************************************************************************
!> \brief Handles CIF (Crystallographic Information File) files
!> \author Teodoro Laino [tlaino]
!> \date 12.2008
! *****************************************************************************
MODULE topology_cif
USE cell_types, ONLY: cell_create,&
cell_release,&
cell_type,&
pbc,&
real_to_scaled,&
scaled_to_real,&
set_cell_param,&
write_cell
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE cp_parser_methods, ONLY: parser_get_next_line,&
parser_get_object,&
parser_search_string
USE cp_parser_types, ONLY: cp_parser_type,&
parser_create,&
parser_release
USE cp_units, ONLY: cp_unit_to_cp2k
USE f77_blas
USE fparser, ONLY: evalf,&
finalizef,&
initf,&
parsef
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: default_string_length,&
dp
USE memory_utilities, ONLY: reallocate
USE string_utilities, ONLY: s2a
USE timings, ONLY: timeset,&
timestop
USE topology_types, ONLY: atom_info_type,&
topology_parameters_type
#include "cp_common_uses.h"
IMPLICIT NONE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_cif'
PRIVATE
PUBLIC :: read_coordinate_cif
CONTAINS
! *****************************************************************************
!> \brief Performs the real task of reading the proper information from the CIF
!> file
!> \par Format Information implemented:
!> _chemical_name
!> _chemical_formula_sum
!> _cell_length_a
!> _cell_length_b
!> _cell_length_c
!> _cell_angle_alpha
!> _cell_angle_beta
!> _cell_angle_gamma
!> _symmetry_space_group_name_h-m
!> _symmetry_equiv_pos_as_xyz
!> _atom_site_label
!> _atom_site_type_symbol
!> _atom_site_fract_x
!> _atom_site_fract_y
!> _atom_site_fract_z
!>
!> \author Teodoro Laino [tlaino]
!> \date 12.2008
! *****************************************************************************
SUBROUTINE read_coordinate_cif (topology,para_env,subsys_section,error)
TYPE(topology_parameters_type) :: topology
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(section_vals_type), POINTER :: subsys_section
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_cif', &
routineP = moduleN//':'//routineN
INTEGER, PARAMETER :: nblock_atom = 10
REAL(KIND=dp), PARAMETER :: realloc_factor = 1.2_dp, &
threshold = 1.0E-6_dp
CHARACTER(LEN=default_string_length) :: s_tag
INTEGER :: handle, ii, iln1, iln2, isym, &
itype, iw, jj, natom, &
natom_orig, newsize
INTEGER, DIMENSION(20) :: ifield
INTEGER, DIMENSION(3) :: periodic
LOGICAL :: check, failure, found, &
my_end, orthorhombic
REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths, r, &
r1, r2, s, s_tmp
TYPE(atom_info_type), POINTER :: atom_info
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_parser_type), POINTER :: parser
NULLIFY(parser, logger)
failure = .FALSE.
logger => cp_error_get_logger(error)
iw = cp_print_key_unit_nr(logger,subsys_section,"PRINT%TOPOLOGY_INFO/CIF_INFO",&
extension=".subsysLog",error=error)
CALL timeset(routineN,handle)
atom_info => topology%atom_info
CALL reallocate(atom_info%label_molname,1,nblock_atom)
CALL reallocate(atom_info%label_resname,1,nblock_atom)
CALL reallocate(atom_info%label_resid,1,nblock_atom)
CALL reallocate(atom_info%label_atmname,1,nblock_atom)
CALL reallocate(atom_info%r,1,3,1,nblock_atom)
CALL reallocate(atom_info%atm_mass,1,nblock_atom)
CALL reallocate(atom_info%atm_charge,1,nblock_atom)
CALL reallocate(atom_info%occup,1,nblock_atom)
CALL reallocate(atom_info%beta,1,nblock_atom)
CALL reallocate(atom_info%element,1,nblock_atom)
IF (iw>0) WRITE(iw,*) " Reading in CIF file ",TRIM(topology%coord_file_name)
CALL parser_create(parser,topology%coord_file_name,&
para_env=para_env,error=error)
! Check for _chemical_name
CALL parser_search_string(parser,"_chemical_name",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
IF (found) THEN
IF (iw>0) WRITE(iw,'(/,A)')" CIF_INFO| _chemical_name :: "//TRIM(parser%input_line(parser%icol:))
END IF
! Check for _chemical_formula_sum
CALL parser_search_string(parser,"_chemical_formula_sum",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
IF (found) THEN
IF (iw>0) WRITE(iw,'(A)')" CIF_INFO| _chemical_formula_sum :: "//TRIM(parser%input_line(parser%icol:))
END IF
! Parsing cell infos
periodic = 1
! Check for _cell_length_a
CALL parser_search_string(parser,"_cell_length_a",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_length_a) was not found in CIF file! "//&
CPSourceFileRef)
CALL cif_get_real(parser,cell_lengths(1),error)
cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1),"angstrom",error=error)
! Check for _cell_length_b
CALL parser_search_string(parser,"_cell_length_b",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_length_b) was not found in CIF file! "//&
CPSourceFileRef)
CALL cif_get_real(parser,cell_lengths(2),error)
cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2),"angstrom",error=error)
! Check for _cell_length_c
CALL parser_search_string(parser,"_cell_length_c",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_length_c) was not found in CIF file! "//&
CPSourceFileRef)
CALL cif_get_real(parser,cell_lengths(3),error)
cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3),"angstrom",error=error)
! Check for _cell_angle_alpha
CALL parser_search_string(parser,"_cell_angle_alpha",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_angle_alpha) was not found in CIF file! "//&
CPSourceFileRef)
CALL cif_get_real(parser,cell_angles(1),error)
cell_angles(1) = cp_unit_to_cp2k(cell_angles(1),"deg",error=error)
! Check for _cell_angle_beta
CALL parser_search_string(parser,"_cell_angle_beta",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_angle_beta) was not found in CIF file! "//&
CPSourceFileRef)
CALL cif_get_real(parser,cell_angles(2),error)
cell_angles(2) = cp_unit_to_cp2k(cell_angles(2),"deg",error=error)
! Check for _cell_angle_gamma
CALL parser_search_string(parser,"_cell_angle_gamma",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_angle_gamma) was not found in CIF file! "//&
CPSourceFileRef)
CALL cif_get_real(parser,cell_angles(3),error)
cell_angles(3) = cp_unit_to_cp2k(cell_angles(3),"deg",error=error)
! Determine if the cell is orthorhombic
orthorhombic = (cell_angles(1)==cp_unit_to_cp2k(90.0_dp,"deg",error=error)).AND.&
(cell_angles(2)==cp_unit_to_cp2k(90.0_dp,"deg",error=error)).AND.&
(cell_angles(3)==cp_unit_to_cp2k(90.0_dp,"deg",error=error))
! Create cell
NULLIFY(cell)
CALL cell_create(cell,periodic=periodic,orthorhombic=orthorhombic,error=error)
CALL set_cell_param(cell,cell_lengths,cell_angles)
CALL write_cell(cell,subsys_section,label="CELL_CIF|",error=error)
! Parse atoms info an fractional coordinates
! Check for _atom_site_label
CALL parser_search_string(parser,"_atom_site_label",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_fatal_level,cp_assertion_failed,routineP,&
"The field (_cell_angle_gamma) was not found in CIF file! "//&
CPSourceFileRef)
ifield = 0
ii = 1
ifield(ii) = 0
CALL parser_get_next_line(parser,1,error=error)
DO WHILE (INDEX(parser%input_line,"_atom_site_")/=0)
ii = ii + 1
CPPostcondition(ii<=20,cp_failure_level,routineP,error,failure)
ifield (ii) = -1
! Check for _atom_site_fract_x
IF (INDEX(parser%input_line,"_atom_site_fract_x")/=0) ifield(ii) = 1
! Check for _atom_site_fract_y
IF (INDEX(parser%input_line,"_atom_site_fract_y")/=0) ifield(ii) = 2
! Check for _atom_site_fract_z
IF (INDEX(parser%input_line,"_atom_site_fract_z")/=0) ifield(ii) = 3
CALL parser_get_next_line(parser,1,error=error)
END DO
! Parse real info
natom = 0
DO WHILE ((INDEX(parser%input_line,"loop_")==0).AND.(parser%input_line(1:1)/="_"))
natom = natom + 1
! Resize in case needed
IF (natom > SIZE(atom_info%label_molname)) THEN
newsize=natom*realloc_factor+nblock_atom
CALL reallocate(atom_info%label_molname,1,newsize)
CALL reallocate(atom_info%label_resname,1,newsize)
CALL reallocate(atom_info%label_resid,1,newsize)
CALL reallocate(atom_info%label_atmname,1,newsize)
CALL reallocate(atom_info%r,1,3,1,newsize)
CALL reallocate(atom_info%atm_mass,1,newsize)
CALL reallocate(atom_info%atm_charge,1,newsize)
CALL reallocate(atom_info%occup,1,newsize)
CALL reallocate(atom_info%beta,1,newsize)
CALL reallocate(atom_info%element,1,newsize)
END IF
ii = 1
itype = 0
DO WHILE (ANY(ifield(ii:)>0))
SELECT CASE (ifield(ii))
CASE (-1)
! Skip this field
CALL parser_get_object (parser,s_tag,error=error)
CASE (0)
CALL parser_get_object (parser,atom_info%label_atmname(natom),error=error)
atom_info%label_molname(natom) = "MOL"//cp_to_string(natom)
atom_info%label_resname(natom) = atom_info%label_molname(natom)
atom_info%label_resid(natom) = 1
atom_info%element(natom) = atom_info%label_atmname(natom)
CASE (1)
CALL cif_get_real(parser,atom_info%r(1,natom),error)
CASE (2)
CALL cif_get_real(parser,atom_info%r(2,natom),error)
CASE (3)
CALL cif_get_real(parser,atom_info%r(3,natom),error)
CASE DEFAULT
! Never reach this point..
CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure)
END SELECT
ii = ii + 1
END DO
s = atom_info%r(1:3,natom)
CALL scaled_to_real(atom_info%r(1:3,natom),s,cell)
CALL parser_get_next_line(parser,1,at_end=my_end,error=error)
IF (my_end) EXIT
END DO
! Preliminary check: check if atoms provided are really unique.. this is a paranoic
! check since they should be REALLY unique.. anyway..
DO ii = 1, natom
r1 = atom_info%r(1:3,ii)
DO jj = ii+1,natom
r2 = atom_info%r(1:3,jj)
r = pbc(r1-r2,cell)
check = (SQRT(DOT_PRODUCT(r,r))>=threshold)
CPPostcondition(check,cp_failure_level,routineP,error,failure)
END DO
END DO
! Parse Symmetry Group and generation elements..
! Check for _symmetry_space_group_name_h-m
CALL parser_search_string(parser,"_symmetry_space_group_name_h-m",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
IF (found) THEN
IF (iw>0) WRITE(iw,'(A)')" CIF_INFO| _symmetry_space_group_name_h-m :: "//TRIM(parser%input_line(parser%icol:))
END IF
! Check for _symmetry_equiv_pos_as_xyz
CALL parser_search_string(parser,"_symmetry_equiv_pos_as_xyz",ignore_case=.FALSE.,found=found,&
begin_line=.FALSE.,search_from_begin_of_file=.TRUE.,error=error)
CALL cp_assert(found,cp_warning_level,cp_assertion_failed,routineP,&
"The field (_symmetry_equiv_pos_as_xyz) was not found in CIF file! "//&
CPSourceFileRef,&
only_ionode=.TRUE.)
IF (iw>0) WRITE(iw,'(A,I0)')" CIF_INFO| Number of atoms before applying symmetry operations :: ",natom
IF (iw>0) WRITE(iw,'(A10,1X,3F12.6)')(TRIM(atom_info%label_atmname(ii)),atom_info%r(1:3,ii),ii=1,natom)
IF (found) THEN
! Apply symmetry elements and generate the whole set of atoms in the unit cell
CALL parser_get_next_line(parser,1,error=error)
isym = 0
natom_orig = natom
DO WHILE ((INDEX(parser%input_line,"loop_")==0).AND.(parser%input_line(1:1)/="_"))
isym = isym + 1
iln1 = INDEX(parser%input_line(1:) ,",")
iln2 = INDEX(parser%input_line(iln1+1:),",")+iln1
CPPostcondition(iln1/=0,cp_failure_level,routineP,error,failure)
CPPostcondition(iln2/=iln1,cp_failure_level,routineP,error,failure)
CALL initf(3)
CALL parsef(1,TRIM(parser%input_line(1:iln1-1)) ,s2a("x","y","z"))
CALL parsef(2,TRIM(parser%input_line(iln1+1:iln2-1)),s2a("x","y","z"))
CALL parsef(3,TRIM(parser%input_line(iln2+1:)) ,s2a("x","y","z"))
Loop_over_unique_atoms: DO ii = 1, natom_orig
CALL real_to_scaled(s_tmp, atom_info%r(1:3,ii), cell)
s(1) = evalf(1,(/s_tmp(1),s_tmp(2),s_tmp(3)/))
s(2) = evalf(2,(/s_tmp(1),s_tmp(2),s_tmp(3)/))
s(3) = evalf(3,(/s_tmp(1),s_tmp(2),s_tmp(3)/))
CALL scaled_to_real(r1, s, cell)
check = .TRUE.
DO jj = 1, natom
r2 = atom_info%r(1:3,jj)
r = pbc(r1-r2,cell)
IF (SQRT(DOT_PRODUCT(r,r))<=threshold) THEN
check = .FALSE.
EXIT
END IF
END DO
! If the atom generated is unique let's add to the atom set..
IF (check) THEN
natom = natom + 1
! Resize in case needed
IF (natom > SIZE(atom_info%label_molname)) THEN
newsize=natom*realloc_factor+nblock_atom
CALL reallocate(atom_info%label_molname,1,newsize)
CALL reallocate(atom_info%label_resname,1,newsize)
CALL reallocate(atom_info%label_resid,1,newsize)
CALL reallocate(atom_info%label_atmname,1,newsize)
CALL reallocate(atom_info%r,1,3,1,newsize)
CALL reallocate(atom_info%atm_mass,1,newsize)
CALL reallocate(atom_info%atm_charge,1,newsize)
CALL reallocate(atom_info%occup,1,newsize)
CALL reallocate(atom_info%beta,1,newsize)
CALL reallocate(atom_info%element,1,newsize)
END IF
atom_info%label_atmname(natom) = atom_info%label_atmname(ii)
atom_info%label_molname(natom) = atom_info%label_molname(ii)
atom_info%label_resname(natom) = atom_info%label_resname(ii)
atom_info%label_resid(natom) = atom_info%label_resid(ii)
atom_info%element(natom) = atom_info%element(ii)
atom_info%r(1:3,natom) = r1
END IF
END DO Loop_over_unique_atoms
CALL finalizef()
CALL parser_get_next_line(parser,1,at_end=my_end,error=error)
IF (my_end) EXIT
END DO
END IF
IF (iw>0) WRITE(iw,'(A,I0)')" CIF_INFO| Number of ssymmetry operations :: ",isym
IF (iw>0) WRITE(iw,'(A,I0)')" CIF_INFO| Number of total atoms :: ",natom
IF (iw>0) WRITE(iw,'(A10,1X,3F12.6)')(TRIM(atom_info%label_atmname(ii)),atom_info%r(1:3,ii),ii=1,natom)
! Releasse local cell type and parser
CALL cell_release(cell,error=error)
CALL parser_release(parser,error=error)
! Reallocate all structures with the exact NATOM size
CALL reallocate(atom_info%label_molname,1,natom)
CALL reallocate(atom_info%label_resname,1,natom)
CALL reallocate(atom_info%label_resid,1,natom)
CALL reallocate(atom_info%label_atmname,1,natom)
CALL reallocate(atom_info%r,1,3,1,natom)
CALL reallocate(atom_info%atm_mass,1,natom)
CALL reallocate(atom_info%atm_charge,1,natom)
CALL reallocate(atom_info%occup,1,natom)
CALL reallocate(atom_info%beta,1,natom)
CALL reallocate(atom_info%element,1,natom)
topology%natoms = natom
topology%molname_generated = .TRUE.
CALL cp_print_key_finished_output(iw,logger,subsys_section,&
"PRINT%TOPOLOGY_INFO/CIF_INFO",error=error)
CALL timestop(handle)
END SUBROUTINE read_coordinate_cif
! *****************************************************************************
!> \brief Reads REAL from the CIF file.. This wrapper is needed in order to
!> treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
!> \author Teodoro Laino [tlaino]
!> \date 12.2008
! *****************************************************************************
SUBROUTINE cif_get_real(parser, r, error)
TYPE(cp_parser_type), POINTER :: parser
REAL(KIND=dp), INTENT(OUT) :: r
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'cif_get_real', &
routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length) :: s_tag
INTEGER :: iln
CALL parser_get_object (parser,s_tag,error=error)
iln = LEN_TRIM(s_tag)
IF (INDEX(s_tag,"(")/=0) iln = INDEX(s_tag,"(")-1
READ(s_tag(1:iln),*) r
END SUBROUTINE cif_get_real
END MODULE topology_cif

6
src/topology_xyz.F
View file

@ -73,7 +73,6 @@ CONTAINS
Frames: DO
! Atom numbers
CALL parser_get_object (parser,natom,error=error)
topology%natoms = natom
frame = frame + 1
IF (frame == 1) THEN
CALL reallocate(atom_info%label_molname, 1, natom)
@ -92,7 +91,6 @@ CONTAINS
END IF
! Dummy line
CALL parser_get_next_line(parser,2,error=error)
topology%molname_generated = .TRUE.
DO j = 1, natom
! Atom coordinates
READ(parser%input_line,*)atom_info%label_atmname(j),&
@ -131,8 +129,8 @@ CPSourceFileRef,&
END DO Frames
CALL parser_release(parser,error=error)
IF(.NOT.topology%para_res) atom_info%label_resid(:) = 1
topology%natoms = natom
topology%natoms = natom
topology%molname_generated = .TRUE.
CALL timestop(handle)
CALL cp_print_key_finished_output(iw,logger,subsys_section,&
"PRINT%TOPOLOGY_INFO/XYZ_INFO",error=error)

61
tests/SE/regtest-3/Periclase.inp Normal file
View file

@ -0,0 +1,61 @@
&GLOBAL
PROJECT MgO
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
&PRINT
&FORCES
&END
&END
METHOD Quickstep
&DFT
&QS
METHOD PM6
&SE
RC_COULOMB [angstrom] 4.175
&END
&END QS
&SCF
MAX_SCF 50
EPS_SCF 1.0E-6
SCF_GUESS mopac
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END
&OUTER_SCF
EPS_SCF 1.0E-6
MAX_SCF 5
&END
&END SCF
&END DFT
&SUBSYS
&CELL
# this is the wrong cell but it's just to test the parse of the CIF
ABC 4.175 4.175 4.175
&END CELL
&TOPOLOGY
COORDINATE CIF
COORD_FILE_NAME ../sample_cif/Periclase.cif
CONNECTIVITY OFF
&DUMP_PDB
&END
&DUMP_PSF
&END
&END
&PRINT
&CELL
UNIT angstrom
&END
&TOPOLOGY_INFO
CIF_INFO
&END
&KINDS
BASIS_SET
POTENTIAL
&END KINDS
&END PRINT
&END SUBSYS
&END FORCE_EVAL

1
tests/SE/regtest-3/TEST_FILES
View file

@ -12,3 +12,4 @@ H2_check.inp 11
O2.inp 11
H2O-3.inp 11
N3-rp_colvar.inp 11
Periclase.inp 11

105
tests/SE/sample_cif/Mercury_phosphide_iodide.cif Executable file
View file

@ -0,0 +1,105 @@
data_1001418
_chemical_name_systematic 'Mercury phosphide iodide (9/5/6)'
_chemical_formula_structural 'Hg9 P5 I6'
_chemical_formula_sum 'Hg9 I6 P5'
_publ_section_title
;
Hg~9~P~5~I~6~: a new mercury(I,II) structural determination
;
loop_
_publ_author_name
'Ledesert, M'
'Rebbah, A'
'Labbe, P'
_journal_name_full
;
Zeitschrift fuer Kristallographie (149,1979-)
;
_journal_coden_astm ZEKRDZ
_journal_volume 192
_journal_year 1990
_journal_page_first 223
_journal_page_last 231
_cell_length_a 13.112(5)
_cell_length_b 12.486(2)
_cell_length_c 17.031(3)
_cell_angle_alpha 90
_cell_angle_beta 119.90(2)
_cell_angle_gamma 90
_cell_volume 2417.1
_cell_formula_units_z 4
_symmetry_space_group_name_h-m 'P 1 21/c 1'
_symmetry_int_tables_number 14
_symmetry_cell_setting monoclinic
loop_
_symmetry_equiv_pos_as_xyz
x,y,z
-x,1/2+y,1/2-z
-x,-y,-z
x,1/2-y,1/2+z
loop_
_atom_type_symbol
_atom_type_oxidation_number
Hg2+ 2.000
Hg1+ 1.000
I1- -1.000
P3- -3.000
P2- -2.000
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_symmetry_multiplicity
_atom_site_Wyckoff_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
_atom_site_attached_hydrogens
_atom_site_calc_flag
Hg1 Hg2+ 2 a 0. 0. 0. 1. 0 d
Hg2 Hg2+ 2 d 0.5 0.5 0. 1. 0 d
Hg3 Hg2+ 4 e 0.0986(3) 0.1764(3) 0.8115(2) 1. 0 d
Hg4 Hg2+ 4 e 0.6452(3) -0.1819(3) 0.7207(2) 1. 0 d
Hg5 Hg2+ 4 e 0.2114(3) -0.2700(4) 0.5343(2) 1. 0 d
Hg6 Hg2+ 4 e 0.4446(3) 0.0579(3) 0.6981(2) 1. 0 d
Hg7 Hg2+ 4 e 0.1250(3) -0.0990(3) 0.7113(2) 1. 0 d
Hg8 Hg2+ 4 e 0.2411(3) -0.3971(3) 0.7642(2) 1. 0 d
Hg9 Hg2+ 4 e 0.1239(4) 0.3140(4) 0.0549(3) 1. 0 d
Hg10 Hg1+ 4 e 0.4005(3) 0.4768(4) 0.4324(3) 1. 0 d
I1 I1- 4 e 0.6122(4) 0.1899(5) 0.2724(4) 1. 0 d
I2 I1- 4 e 0.3115(5) 0.2860(5) -0.0247(4) 1. 0 d
I3 I1- 4 e 0.3449(5) -0.2509(6) 0.9575(4) 1. 0 d
I4 I1- 4 e -0.0471(4) -0.1122(5) 0.8060(3) 1. 0 d
I5 I1- 4 e 0.0936(6) -0.0285(6) 0.4404(5) 1. 0 d
I6 I1- 4 e 0.7068(4) 0.5153(5) 0.4114(4) 1. 0 d
P1 P3- 4 e 0.584(1) 0.505(2) 0.160(1) 1. 0 d
P2 P2- 4 e 0.928(2) 0.289(2) 0.325(1) 1. 0 d
P3 P2- 4 e 0.291(2) 0.355(2) 0.208(1) 1. 0 d
P4 P2- 4 e 0.965(2) 0.307(2) 0.459(1) 1. 0 d
P5 P2- 4 e 0.226(1) 0.437(2) 0.286(1) 1. 0 d
loop_
_atom_site_aniso_label
_atom_site_aniso_U_11
_atom_site_aniso_U_12
_atom_site_aniso_U_13
_atom_site_aniso_U_22
_atom_site_aniso_U_23
_atom_site_aniso_U_33
Hg1 0.050(2) 0.011(3) 0.027(2) 0.014(3) 0.010(3) 0.034(2)
Hg2 0.028(2) 0.010(3) 0.010(1) 0.041(4) 0.006(3) 0.019(2)
Hg3 0.032(1) 0.012(2) 0.0319(9) 0.032(2) 0.005(2) 0.046(1)
Hg4 0.039(1) 0.006(2) 0.030(1) 0.015(2) 0.007(2) 0.058(2)
Hg5 0.043(1) -0.014(2) 0.0447(9) 0.052(3) -0.007(2) 0.062(2)
Hg6 0.030(1) 0.005(2) 0.0331(9) 0.025(2) -0.004(2) 0.055(2)
Hg7 0.049(2) -0.013(2) 0.032(1) 0.021(2) -0.016(2) 0.043(2)
Hg8 0.022(1) 0.004(2) 0.015(1) 0.031(2) 0.005(2) 0.039(2)
Hg9 0.060(2) -0.043(3) 0.021(2) 0.083(3) -0.037(3) 0.044(2)
Hg10 0.024(2) -0.001(2) 0.007(1) 0.058(3) -0.026(2) 0.042(2)
I1 0.017(2) .000(3) 0.019(1) 0.018(3) .000(3) 0.043(3)
I2 0.033(2) -0.010(3) 0.026(2) 0.031(4) -0.002(3) 0.048(3)
I3 0.053(2) -0.001(3) 0.046(2) 0.036(4) -0.006(3) 0.058(3)
I4 0.033(2) 0.005(3) 0.021(1) 0.018(3) 0.001(3) 0.026(2)
I5 0.049(4) 0.005(4) 0.014(3) 0.032(5) 0.015(4) 0.076(5)
I6 0.043(2) 0.013(3) 0.034(2) 0.029(4) 0.008(3) 0.049(3)
_refine_ls_r_factor_all 0.046

226
tests/SE/sample_cif/Periclase.cif Executable file
View file

@ -0,0 +1,226 @@
data_9006799
_chemical_name Periclase
loop_
_publ_author_name
'Zhang J'
_journal_name_full 'Physics and Chemistry of Minerals'
_journal_volume 27
_journal_year 2000
_journal_page_first 145
_journal_page_last 148
_publ_section_title
;
Effect of pressure on the thermal expansion of MgO up to 8.2 GPa
Sample: Run2 at T = 473 K, P = 5.44 GPa
;
_chemical_formula_sum 'Mg O'
_cell_length_a 4.1749
_cell_length_b 4.1749
_cell_length_c 4.1749
_cell_angle_alpha 90
_cell_angle_beta 90
_cell_angle_gamma 90
_cell_volume 72.768
_symmetry_space_group_name_h-m 'F m 3 m'
loop_
_symmetry_equiv_pos_as_xyz
x,y,z
x,1/2+y,1/2+z
1/2+x,y,1/2+z
1/2+x,1/2+y,z
z,-x,y
z,1/2-x,1/2+y
1/2+z,-x,1/2+y
1/2+z,1/2-x,y
-y,z,-x
-y,1/2+z,1/2-x
1/2-y,z,1/2-x
1/2-y,1/2+z,-x
x,-y,z
x,1/2-y,1/2+z
1/2+x,-y,1/2+z
1/2+x,1/2-y,z
-z,x,-y
-z,1/2+x,1/2-y
1/2-z,x,1/2-y
1/2-z,1/2+x,-y
y,-z,x
y,1/2-z,1/2+x
1/2+y,-z,1/2+x
1/2+y,1/2-z,x
-x,y,-z
-x,1/2+y,1/2-z
1/2-x,y,1/2-z
1/2-x,1/2+y,-z
x,-z,-y
x,1/2-z,1/2-y
1/2+x,-z,1/2-y
1/2+x,1/2-z,-y
-z,y,x
-z,1/2+y,1/2+x
1/2-z,y,1/2+x
1/2-z,1/2+y,x
y,-x,-z
y,1/2-x,1/2-z
1/2+y,-x,1/2-z
1/2+y,1/2-x,-z
-x,z,y
-x,1/2+z,1/2+y
1/2-x,z,1/2+y
1/2-x,1/2+z,y
z,-y,-x
z,1/2-y,1/2-x
1/2+z,-y,1/2-x
1/2+z,1/2-y,-x
-y,x,z
-y,1/2+x,1/2+z
1/2-y,x,1/2+z
1/2-y,1/2+x,z
x,z,y
x,1/2+z,1/2+y
1/2+x,z,1/2+y
1/2+x,1/2+z,y
-z,-y,-x
-z,1/2-y,1/2-x
1/2-z,-y,1/2-x
1/2-z,1/2-y,-x
y,x,z
y,1/2+x,1/2+z
1/2+y,x,1/2+z
1/2+y,1/2+x,z
-x,-z,-y
-x,1/2-z,1/2-y
1/2-x,-z,1/2-y
1/2-x,1/2-z,-y
z,y,x
z,1/2+y,1/2+x
1/2+z,y,1/2+x
1/2+z,1/2+y,x
-y,-x,-z
-y,1/2-x,1/2-z
1/2-y,-x,1/2-z
1/2-y,1/2-x,-z
z,x,-y
z,1/2+x,1/2-y
1/2+z,x,1/2-y
1/2+z,1/2+x,-y
-y,-z,x
-y,1/2-z,1/2+x
1/2-y,-z,1/2+x
1/2-y,1/2-z,x
x,y,-z
x,1/2+y,1/2-z
1/2+x,y,1/2-z
1/2+x,1/2+y,-z
-z,-x,y
-z,1/2-x,1/2+y
1/2-z,-x,1/2+y
1/2-z,1/2-x,y
y,z,-x
y,1/2+z,1/2-x
1/2+y,z,1/2-x
1/2+y,1/2+z,-x
-x,-y,z
-x,1/2-y,1/2+z
1/2-x,-y,1/2+z
1/2-x,1/2-y,z
-z,x,y
-z,1/2+x,1/2+y
1/2-z,x,1/2+y
1/2-z,1/2+x,y
y,-z,-x
y,1/2-z,1/2-x
1/2+y,-z,1/2-x
1/2+y,1/2-z,-x
-x,y,z
-x,1/2+y,1/2+z
1/2-x,y,1/2+z
1/2-x,1/2+y,z
z,-x,-y
z,1/2-x,1/2-y
1/2+z,-x,1/2-y
1/2+z,1/2-x,-y
-y,z,x
-y,1/2+z,1/2+x
1/2-y,z,1/2+x
1/2-y,1/2+z,x
x,-y,-z
x,1/2-y,1/2-z
1/2+x,-y,1/2-z
1/2+x,1/2-y,-z
-x,z,-y
-x,1/2+z,1/2-y
1/2-x,z,1/2-y
1/2-x,1/2+z,-y
z,-y,x
z,1/2-y,1/2+x
1/2+z,-y,1/2+x
1/2+z,1/2-y,x
-y,x,-z
-y,1/2+x,1/2-z
1/2-y,x,1/2-z
1/2-y,1/2+x,-z
x,-z,y
x,1/2-z,1/2+y
1/2+x,-z,1/2+y
1/2+x,1/2-z,y
-z,y,-x
-z,1/2+y,1/2-x
1/2-z,y,1/2-x
1/2-z,1/2+y,-x
y,-x,z
y,1/2-x,1/2+z
1/2+y,-x,1/2+z
1/2+y,1/2-x,z
-x,-z,y
-x,1/2-z,1/2+y
1/2-x,-z,1/2+y
1/2-x,1/2-z,y
z,y,-x
z,1/2+y,1/2-x
1/2+z,y,1/2-x
1/2+z,1/2+y,-x
-y,-x,z
-y,1/2-x,1/2+z
1/2-y,-x,1/2+z
1/2-y,1/2-x,z
x,z,-y
x,1/2+z,1/2-y
1/2+x,z,1/2-y
1/2+x,1/2+z,-y
-z,-y,x
-z,1/2-y,1/2+x
1/2-z,-y,1/2+x
1/2-z,1/2-y,x
y,x,-z
y,1/2+x,1/2-z
1/2+y,x,1/2-z
1/2+y,1/2+x,-z
-z,-x,-y
-z,1/2-x,1/2-y
1/2-z,-x,1/2-y
1/2-z,1/2-x,-y
y,z,x
y,1/2+z,1/2+x
1/2+y,z,1/2+x
1/2+y,1/2+z,x
-x,-y,-z
-x,1/2-y,1/2-z
1/2-x,-y,1/2-z
1/2-x,1/2-y,-z
z,x,y
z,1/2+x,1/2+y
1/2+z,x,1/2+y
1/2+z,1/2+x,y
-y,-z,-x
-y,1/2-z,1/2-x
1/2-y,-z,1/2-x
1/2-y,1/2-z,-x
loop_
_atom_site_label
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
Mg 0.00000 0.00000 0.00000
O 0.50000 0.50000 0.50000