From 39d24ab3ce2a2d3a74842b50e799d47a8d78a26b Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Wed, 31 Dec 2008 08:38:47 +0000 Subject: [PATCH] Added CIF (Christallographic Information File) as possible format for inputting 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 --- src/OBJECTDEFS | 1 + src/input_constants.F | 11 +- src/input_cp2k_subsys.F | 38 +- src/topology.F | 31 +- src/topology_cif.F | 419 ++++++++++++++++++ src/topology_xyz.F | 6 +- tests/SE/regtest-3/Periclase.inp | 61 +++ tests/SE/regtest-3/TEST_FILES | 1 + .../sample_cif/Mercury_phosphide_iodide.cif | 105 +++++ tests/SE/sample_cif/Periclase.cif | 226 ++++++++++ 10 files changed, 861 insertions(+), 38 deletions(-) create mode 100644 src/topology_cif.F create mode 100644 tests/SE/regtest-3/Periclase.inp create mode 100755 tests/SE/sample_cif/Mercury_phosphide_iodide.cif create mode 100755 tests/SE/sample_cif/Periclase.cif diff --git a/src/OBJECTDEFS b/src/OBJECTDEFS index 9a3ba2854c..c39677d596 100644 --- a/src/OBJECTDEFS +++ b/src/OBJECTDEFS @@ -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\ diff --git a/src/input_constants.F b/src/input_constants.F index 6e7af42bdb..37765ad8c4 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -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,& diff --git a/src/input_cp2k_subsys.F b/src/input_cp2k_subsys.F index 1b412d8a4a..89e53690c9 100644 --- a/src/input_cp2k_subsys.F +++ b/src/input_cp2k_subsys.F @@ -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) diff --git a/src/topology.F b/src/topology.F index eabd4dd36f..a4b5072e8b 100644 --- a/src/topology.F +++ b/src/topology.F @@ -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) diff --git a/src/topology_cif.F b/src/topology_cif.F new file mode 100644 index 0000000000..cd85fbb4e3 --- /dev/null +++ b/src/topology_cif.F @@ -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 diff --git a/src/topology_xyz.F b/src/topology_xyz.F index a4a5a47afb..766f3cb0c4 100644 --- a/src/topology_xyz.F +++ b/src/topology_xyz.F @@ -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) diff --git a/tests/SE/regtest-3/Periclase.inp b/tests/SE/regtest-3/Periclase.inp new file mode 100644 index 0000000000..05d2e59c02 --- /dev/null +++ b/tests/SE/regtest-3/Periclase.inp @@ -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 diff --git a/tests/SE/regtest-3/TEST_FILES b/tests/SE/regtest-3/TEST_FILES index 2c9def6c19..7e866a9602 100644 --- a/tests/SE/regtest-3/TEST_FILES +++ b/tests/SE/regtest-3/TEST_FILES @@ -12,3 +12,4 @@ H2_check.inp 11 O2.inp 11 H2O-3.inp 11 N3-rp_colvar.inp 11 +Periclase.inp 11 diff --git a/tests/SE/sample_cif/Mercury_phosphide_iodide.cif b/tests/SE/sample_cif/Mercury_phosphide_iodide.cif new file mode 100755 index 0000000000..d206baf9e3 --- /dev/null +++ b/tests/SE/sample_cif/Mercury_phosphide_iodide.cif @@ -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 diff --git a/tests/SE/sample_cif/Periclase.cif b/tests/SE/sample_cif/Periclase.cif new file mode 100755 index 0000000000..a12be7f993 --- /dev/null +++ b/tests/SE/sample_cif/Periclase.cif @@ -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