TABLE OF CONTENTS
- ABINIT/m_effective_potential_file
- m_effective_potential_file/coeffs_xml2effpot
- m_effective_potential_file/effective_potential_file_getDimCoeff
- m_effective_potential_file/effective_potential_file_getDimMD
- m_effective_potential_file/effective_potential_file_getDimStrainCoupling
- m_effective_potential_file/effective_potential_file_getDimSystem
- m_effective_potential_file/effective_potential_file_getType
- m_effective_potential_file/effective_potential_file_mapHistToRef
- m_effective_potential_file/effective_potential_file_readMDfile
- m_effective_potential_file/elementfromline
- m_effective_potential_file/rdfromline
- m_effective_potential_file/rdfromline_value
- m_effective_potential_file/rmtabfromline
- m_effective_potential_file/system_ddb2effpot
- m_effective_potential_file/system_getDimFromXML
- m_effective_potential_file/system_xml2effpot
- m_effpot_xml/char_c2f
- m_effpot_xml/char_f2c
ABINIT/m_effective_potential_file [ Modules ]
NAME
m_effective_potential_file
FUNCTION
This module contains all routine to read the effective potential from files Can also read coefficients from XML (XML or DDB)
COPYRIGHT
Copyright (C) 2000-2024 ABINIT group (AM) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_effective_potential_file 26 27 use defs_basis 28 use m_errors 29 use m_abicore 30 use m_xmpi 31 use m_harmonics_terms 32 use m_anharmonics_terms 33 use m_effective_potential 34 use m_ifc 35 use m_ddb 36 use m_ddb_hdr 37 #if defined HAVE_NETCDF 38 use netcdf 39 #endif 40 41 use m_io_tools, only : open_file, get_unit 42 use m_geometry, only : xcart2xred, xred2xcart, metric 43 use m_symfind, only : symfind, symlatt 44 use m_crystal, only : crystal_t, crystal_init 45 use m_dynmat, only : dfpt_prtph 46 use m_abihist, only : abihist,abihist_init,abihist_free,abihist_copy,read_md_hist 47 use m_ddb_internalstr, only : ddb_internalstr 48 49 implicit none 50 51 public :: effective_potential_file_getDimCoeff 52 public :: effective_potential_file_getDimMD 53 public :: effective_potential_file_getDimSystem 54 public :: effective_potential_file_getDimStrainCoupling 55 public :: effective_potential_file_getType 56 public :: effective_potential_file_mapHistToRef 57 public :: effective_potential_file_read 58 public :: effective_potential_file_readDisplacement 59 public :: effective_potential_file_readMDfile 60 private :: coeffs_xml2effpot 61 private :: system_getDimFromXML 62 private :: system_xml2effpot 63 private :: system_ddb2effpot 64 65 #ifndef HAVE_XML 66 private :: rdfromline 67 private :: rmtabfromline 68 private :: rdfromline_value 69 private :: elementfromline 70 #endif 71 72 #if defined HAVE_XML 73 public :: effpot_xml_checkXML 74 public :: effpot_xml_getDimCoeff 75 public :: effpot_xml_readSystem 76 public :: effpot_xml_getValue 77 public :: effpot_xml_getAttribute 78 public :: effpot_xml_getDimSystem 79 80 interface 81 subroutine effpot_xml_readSystem(filename,natom,& 82 & ntypat,nrpt,nqpt,amu,atmfrc,cell,dynmat,elastic_constants,& 83 & energy,epsilon_inf,ewald_atmfrc,& 84 & phfrq,rprimd,qph1l,short_atmfrc,typat,xcart,zeff)& 85 & bind(C,name="effpot_xml_readSystem") 86 use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT 87 integer(C_INT) :: natom,ntypat,nrpt,nqpt 88 integer(C_INT) :: typat(natom) 89 integer(C_INT) :: cell(3,nrpt) 90 real(C_DOUBLE) :: energy 91 real(C_DOUBLE) :: dynmat(2,3,natom,3,natom,nqpt) 92 real(C_DOUBLE) :: phfrq(3*natom,nqpt),qph1l(3,nqpt) 93 real(C_DOUBLE) :: atmfrc(3,natom,3,natom,nrpt) 94 real(C_DOUBLE) :: short_atmfrc(3,natom,3,natom,nrpt) 95 real(C_DOUBLE) :: ewald_atmfrc(3,natom,3,natom,nrpt) 96 real(C_DOUBLE) :: amu(ntypat),rprimd(3,3),epsilon_inf(3,3) 97 real(C_DOUBLE) :: zeff(3,3,natom) 98 real(C_DOUBLE) :: elastic_constants(6,6),xcart(3,natom) 99 character(kind=C_CHAR) :: filename(*) 100 end subroutine effpot_xml_readSystem 101 end interface 102 103 interface 104 subroutine effpot_xml_readStrainCoupling(filename,natom,& 105 & nrpt,voigt,elastic3rd,elastic_displacement,& 106 & strain_coupling,phonon_strain_atmfrc,phonon_straincell)& 107 & bind(C,name="effpot_xml_readStrainCoupling") 108 use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT 109 integer(C_INT) :: natom 110 integer(C_INT) :: nrpt,voigt 111 integer(c_INT) :: phonon_straincell(3,nrpt) 112 real(C_DOUBLE) :: elastic3rd(6,6),elastic_displacement(6,3,natom) 113 real(C_DOUBLE) :: strain_coupling(3,natom) 114 real(C_DOUBLE) :: phonon_strain_atmfrc(3,natom,3,natom,nrpt) 115 character(kind=C_CHAR) :: filename(*) 116 end subroutine effpot_xml_readStrainCoupling 117 end interface 118 119 interface 120 subroutine effpot_xml_readCoeff(filename,ncoeff,ndisp,nterm,& 121 & coefficient,atindx,cell,direction,power_disp,& 122 & power_strain,strain,weight)& 123 & bind(C,name="effpot_xml_readCoeff") 124 use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT 125 character(kind=C_CHAR) :: filename(*) 126 integer(C_INT) :: ncoeff,ndisp,nterm 127 integer(C_INT) :: atindx(ncoeff,nterm,2,ndisp) 128 integer(C_INT) :: cell(ncoeff,nterm,3,2,ndisp) 129 integer(C_INT) :: direction(ncoeff,nterm,ndisp) 130 integer(C_INT) :: strain(ncoeff,nterm,ndisp) 131 integer(C_INT) :: power_disp(ncoeff,nterm,ndisp) 132 integer(C_INT) :: power_strain(ncoeff,nterm,ndisp) 133 real(C_DOUBLE) :: coefficient(ncoeff) 134 real(C_DOUBLE) :: weight(ncoeff,nterm) 135 end subroutine effpot_xml_readCoeff 136 end interface 137 138 interface 139 subroutine effpot_xml_getDimSystem(filename,natom,ntypat,nqpt,nrpt1,nrpt2)& 140 & bind(C,name="effpot_xml_getDimSystem") 141 use, intrinsic :: iso_c_binding, only : C_CHAR,C_INT 142 integer(C_INT) :: natom,ntypat,nqpt,nrpt1,nrpt2 143 character(kind=C_CHAR) :: filename(*) 144 end subroutine effpot_xml_getDimSystem 145 end interface 146 147 interface 148 subroutine effpot_xml_getDimStrainCoupling(filename,nrpt,voigt)& 149 & bind(C,name="effpot_xml_getDimStrainCoupling") 150 use, intrinsic :: iso_c_binding, only : C_CHAR,C_INT 151 integer(C_INT) :: voigt 152 integer(C_INT) :: nrpt 153 character(kind=C_CHAR) :: filename(*) 154 end subroutine effpot_xml_getDimStrainCoupling 155 end interface 156 157 interface 158 subroutine effpot_xml_getDimCoeff(filename,ncoeff,nterm_max,ndisp_max)& 159 & bind(C,name="effpot_xml_getDimCoeff") 160 use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT,C_PTR 161 character(kind=C_CHAR) :: filename(*) 162 ! character(kind=C_CHAR) :: name(*) 163 type(C_PTR) :: name 164 integer(C_INT) :: ncoeff,ndisp_max,nterm_max 165 end subroutine effpot_xml_getDimCoeff 166 end interface 167 168 169 interface 170 subroutine effpot_xml_checkXML(filename,name_root) & 171 & bind(C,name="effpot_xml_checkXML") 172 use, intrinsic :: iso_c_binding, only : C_CHAR 173 character(kind=C_CHAR) :: filename(*),name_root(*) 174 end subroutine effpot_xml_checkXML 175 end interface 176 177 interface 178 subroutine effpot_xml_getValue(filename,name_value,value_result) & 179 & bind(C,name="effpot_xml_getValue") 180 use, intrinsic :: iso_c_binding, only : C_CHAR 181 implicit none 182 character(kind=C_CHAR) :: filename(*),name_value(*) 183 character(kind=C_CHAR) :: value_result 184 end subroutine effpot_xml_getValue 185 end interface 186 187 interface 188 subroutine effpot_xml_getAttribute(filename,name_value,name_attribute) & 189 & bind(C,name="effpot_xml_getAttribute") 190 use, intrinsic :: iso_c_binding, only : C_CHAR 191 character(kind=C_CHAR) :: filename(*),name_value(*),name_attribute(*) 192 end subroutine effpot_xml_getAttribute 193 end interface 194 195 interface 196 subroutine effpot_xml_getNumberKey(filename,name_value,number) & 197 & bind(C,name="effpot_xml_getNumberKey") 198 use, intrinsic :: iso_c_binding, only : C_CHAR,C_INT 199 character(kind=C_CHAR) :: filename(*),name_value(*) 200 integer(C_INT) :: number 201 end subroutine effpot_xml_getNumberKey 202 end interface 203 204 #endif 205 206 integer,parameter :: XML_RECL = 50000
m_effective_potential_file/coeffs_xml2effpot [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
coeffs_xml2effpot
FUNCTION
Open xml file of effective potentiel, then reads the variables and store them in effective potentential type
INPUTS
filename = path of input or output file comm=MPI communicator
OUTPUT
eff_pot<type(effective_potential_type)> = effective_potential datatype
SOURCE
2837 subroutine coeffs_xml2effpot(eff_pot,filename,comm) 2838 2839 use m_atomdata 2840 use m_polynomial_coeff 2841 use m_polynomial_term 2842 use m_crystal, only : symbols_crystal 2843 #if defined HAVE_XML 2844 use, intrinsic :: iso_c_binding, only : C_CHAR,C_PTR,c_f_pointer 2845 #endif 2846 2847 !Arguments ------------------------------------ 2848 !scalars 2849 character(len=*),intent(in) :: filename 2850 integer, intent(in) :: comm 2851 !arrays 2852 type(effective_potential_type), intent(inout) :: eff_pot 2853 2854 !Local variables------------------------------- 2855 !scalar 2856 integer :: ii,jj,my_rank,ndisp,ncoeff,nterm_max,nstrain,ndisp_max,nproc,nterm 2857 ! character(len=200),allocatable :: name(:) 2858 character(len=200) :: name 2859 #ifdef HAVE_XML 2860 integer :: icoeff,iterm 2861 #endif 2862 2863 #ifndef HAVE_XML 2864 integer :: funit = 1,ios = 0 2865 integer :: icoeff,idisp,istrain,iterm,mu 2866 logical :: found,found2,displacement 2867 character (len=XML_RECL) :: line,readline 2868 character (len=XML_RECL) :: strg,strg1 2869 #endif 2870 character(len=500) :: message 2871 character(len=264) :: filename_tmp 2872 character(len=5),allocatable :: symbols(:) 2873 integer,parameter :: master=0 2874 logical :: iam_master 2875 logical :: debug 2876 !arrays 2877 real(dp),allocatable :: coefficient(:),weight(:,:) 2878 integer,allocatable :: atindx(:,:,:,:), cell(:,:,:,:,:),direction(:,:,:),power_disp(:,:,:) 2879 integer,allocatable :: strain(:,:,:),power_strain(:,:,:) 2880 type(polynomial_coeff_type),dimension(:),allocatable :: coeffs 2881 type(polynomial_term_type),dimension(:,:),allocatable :: terms 2882 2883 ! ************************************************************************* 2884 2885 filename_tmp = trim(filename) 2886 !Open the atomicdata XML file for reading 2887 write(message,'(a,a)')'-Opening the file ',filename_tmp 2888 2889 call wrtout(ab_out,message,'COLL') 2890 call wrtout(std_out,message,'COLL') 2891 2892 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2893 iam_master = (my_rank == master) 2894 2895 !Get Dimention of system and allocation/initialisation of array 2896 ncoeff = 0 2897 nterm = 0 2898 ndisp = 0 2899 nstrain = 0 2900 call effective_potential_file_getDimCoeff(filename,ncoeff,ndisp_max,nterm_max) 2901 2902 ! Do some checks 2903 if (nterm_max<=0) then 2904 write(message, '(a,a,a)' )& 2905 & ' Unable to read the number of terms in ',trim(filename),ch10 2906 ABI_ERROR(message) 2907 end if 2908 2909 if (ndisp_max<=0) then 2910 write(message, '(a,a,a)' )& 2911 & ' Unable to read the number of displacement in ',trim(filename),ch10 2912 ABI_ERROR(message) 2913 end if 2914 2915 !Allocation ov the polynomial coeff type 2916 ABI_MALLOC(coeffs,(ncoeff)) 2917 2918 if(iam_master)then 2919 2920 #if defined HAVE_XML 2921 write(message,'(3a)')'-Reading the file ',trim(filename),& 2922 & ' with LibXML library' 2923 #else 2924 write(message,'(3a)')'-Reading the file ',trim(filename),& 2925 & ' with Fortran' 2926 #endif 2927 call wrtout(ab_out,message,'COLL') 2928 call wrtout(std_out,message,'COLL') 2929 2930 2931 ABI_MALLOC(symbols,(eff_pot%crystal%natom)) 2932 ! Get the symbols arrays 2933 call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,& 2934 & eff_pot%crystal%npsp,symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl) 2935 2936 2937 !Read with libxml librarie 2938 #if defined HAVE_XML 2939 2940 ABI_MALLOC(terms,(ncoeff,nterm_max)) 2941 ABI_MALLOC(atindx,(ncoeff,nterm_max,2,ndisp_max)) 2942 ABI_MALLOC(coefficient,(ncoeff)) 2943 ABI_MALLOC(cell,(ncoeff,nterm_max,3,2,ndisp_max)) 2944 ABI_MALLOC(direction,(ncoeff,nterm_max,ndisp_max)) 2945 ABI_MALLOC(strain,(ncoeff,nterm_max,ndisp_max)) 2946 ABI_MALLOC(power_disp,(ncoeff,nterm_max,ndisp_max)) 2947 ABI_MALLOC(power_strain,(ncoeff,nterm_max,ndisp_max)) 2948 ABI_MALLOC(weight,(ncoeff,nterm_max)) 2949 2950 ! Read the values of this term with libxml 2951 call effpot_xml_readCoeff(char_f2c(trim(filename)),ncoeff,ndisp_max,nterm_max,& 2952 & coefficient,atindx,cell,direction,power_disp,power_strain,& 2953 & strain,weight) 2954 ! In the XML the atom index begin to zero 2955 ! Need to shift for fortran array 2956 atindx(:,:,:,:) = atindx(:,:,:,:) + 1 2957 2958 do icoeff=1,ncoeff 2959 do iterm=1,nterm_max 2960 ! Initialisation of the polynomial_term structure with the values from the 2961 call polynomial_term_init(atindx(icoeff,iterm,:,:),cell(icoeff,iterm,:,:,:),& 2962 & direction(icoeff,iterm,:),ndisp_max,ndisp_max,terms(icoeff,iterm),& 2963 & power_disp(icoeff,iterm,:),power_strain(icoeff,iterm,:),& 2964 & strain(icoeff,iterm,:),weight(icoeff,iterm),check=.true.) 2965 end do 2966 ! Initialisation of the polynomial_coefficent structure with the values 2967 call polynomial_coeff_init(coefficient(icoeff),nterm_max,coeffs(icoeff),& 2968 & terms(icoeff,:),check=.true.) 2969 2970 ! Get the name of this coefficient and set it 2971 ! Try to find the index of the term corresponding to the interation in the 2972 ! reference cell (000) in order to compute the name correctly... 2973 ! If this coeff is not in the ref cell, take by default the first term 2974 if(coeffs(icoeff)%nterm > 0)then 2975 call polynomial_coeff_getName(name,coeffs(icoeff),symbols,recompute=.true.) 2976 call polynomial_coeff_setName(name,coeffs(icoeff)) 2977 end if 2978 2979 ! Free them all 2980 do iterm=1,nterm_max 2981 call polynomial_term_free(terms(icoeff,iterm)) 2982 end do 2983 end do 2984 2985 #else 2986 ABI_MALLOC(terms,(1,nterm_max)) 2987 ABI_MALLOC(atindx,(1,1,2,ndisp_max)) 2988 ABI_MALLOC(coefficient,(1)) 2989 ABI_MALLOC(cell,(1,1,3,2,ndisp_max)) 2990 ABI_MALLOC(direction,(1,1,ndisp_max)) 2991 ABI_MALLOC(strain,(1,1,ndisp_max)) 2992 ABI_MALLOC(power_disp,(1,1,ndisp_max)) 2993 ABI_MALLOC(power_strain,(1,1,ndisp_max)) 2994 ABI_MALLOC(weight,(1,1)) 2995 ! Loop over the file 2996 ! Read the values of all the terms with fortran 2997 if (open_file(filename,message,unit=funit,form="formatted",& 2998 & status="old",action="read") /= 0) then 2999 ABI_ERROR(message) 3000 end if 3001 3002 ! Start a reading loop in fortran 3003 rewind(unit=funit) 3004 ios = 0 3005 found=.false. 3006 3007 ! Initialisation of counter 3008 icoeff = 0 3009 3010 ! Parser 3011 do while (ios==0) 3012 read(funit,'(a)',iostat=ios) readline 3013 if (ios == 0) then 3014 call rmtabfromline(readline) 3015 line=adjustl(readline) 3016 if ((line(1:12)==char(60)//'coefficient')) then 3017 ! Read headers of coefficient 3018 call rdfromline('text',line,strg) 3019 if (strg/="") then 3020 name=trim(strg) 3021 end if 3022 call rdfromline('value',line,strg) 3023 if (strg/="") then 3024 strg1=trim(strg) 3025 read(strg1,*) coefficient(1) 3026 else 3027 coefficient(1) = zero 3028 end if 3029 ! End read headers of coefficient 3030 ! Reset counter 3031 found = .false. 3032 atindx = 0; cell = 0 ; direction = 0 3033 strain = 0; power_strain = 0; power_disp = 0 3034 iterm = 0 3035 idisp = 0 3036 istrain = 0 3037 nterm = 0 3038 do while (.not.found) 3039 read(funit,'(a)',iostat=ios) readline 3040 call rmtabfromline(readline) 3041 line=adjustl(readline) 3042 if ((line(1:13)==char(60)//'/coefficient')) then 3043 found= .true. 3044 cycle 3045 end if 3046 if ((line(1:5)==char(60)//'term')) then 3047 nterm = nterm + 1 3048 ndisp = 0 3049 nstrain = 0 3050 idisp = 0 3051 istrain = 0 3052 displacement = .true. 3053 call rdfromline('weight',line,strg) 3054 if (strg/="") then 3055 strg1=trim(strg) 3056 read(strg1,*) weight 3057 end if 3058 do while(displacement) 3059 read(funit,'(a)',iostat=ios) readline 3060 call rmtabfromline(readline) 3061 line=adjustl(readline) 3062 if ((line(1:6)==char(60)//'/term')) then 3063 displacement = .false. 3064 end if 3065 if ((line(1:7)==char(60)//'strain')) then 3066 nstrain = nstrain + 1 3067 istrain = istrain + 1 3068 call rdfromline('power',line,strg) 3069 if (strg/="") then 3070 strg1=trim(strg) 3071 read(strg1,*) power_strain(1,1,istrain) 3072 end if 3073 call rdfromline('voigt',line,strg) 3074 if (strg/="") then 3075 strg1=trim(strg) 3076 read(strg1,*) strain(1,1,istrain) 3077 end if 3078 end if 3079 if ((line(1:18)==char(60)//'displacement_diff')) then 3080 ndisp = ndisp + 1 3081 idisp = idisp + 1 3082 found2=.true. 3083 call rdfromline('atom_a',line,strg) 3084 if (strg/="") then 3085 strg1=trim(strg) 3086 read(strg1,*) atindx(1,1,1,idisp) 3087 end if 3088 call rdfromline('atom_b',line,strg) 3089 if (strg/="") then 3090 strg1=trim(strg) 3091 read(strg1,*) atindx(1,1,2,idisp) 3092 end if 3093 call rdfromline('direction',line,strg) 3094 if (strg/="") then 3095 strg1=trim(strg) 3096 if (trim(strg1).eq."x") direction(1,1,idisp) = 1 3097 if (trim(strg1).eq."y") direction(1,1,idisp) = 2 3098 if (trim(strg1).eq."z") direction(1,1,idisp) = 3 3099 end if 3100 call rdfromline('power',line,strg) 3101 if (strg/="") then 3102 strg1=trim(strg) 3103 read(strg1,*) power_disp(1,1,idisp) 3104 end if 3105 do while(found2) 3106 read(funit,'(a)',iostat=ios) readline 3107 call rmtabfromline(readline) 3108 line=adjustl(readline) 3109 if ((line(1:7)==char(60)//'cell_a')) then 3110 call rdfromline_value('cell_a',line,strg) 3111 if (strg/="") then 3112 strg1=trim(strg) 3113 read(strg1,*) (cell(1,1,mu,1,idisp),mu=1,3) 3114 else 3115 read(funit,'(a)',iostat=ios) readline 3116 call rmtabfromline(readline) 3117 line=adjustl(readline) 3118 call rdfromline_value('cell_a',line,strg) 3119 if (strg/="") then 3120 strg1=trim(strg) 3121 read(strg1,*)(cell(1,1,mu,1,idisp),mu=1,3) 3122 else 3123 strg1=trim(line) 3124 read(strg1,*)(cell(1,1,mu,1,idisp),mu=1,3) 3125 end if 3126 end if 3127 end if 3128 if ((line(1:7)==char(60)//'cell_b')) then 3129 call rdfromline_value('cell_b',line,strg) 3130 if (strg/="") then 3131 strg1=trim(strg) 3132 read(strg1,*) (cell(1,1,mu,2,idisp),mu=1,3) 3133 else 3134 read(funit,'(a)',iostat=ios) readline 3135 call rmtabfromline(readline) 3136 line=adjustl(readline) 3137 call rdfromline_value('cell_b',line,strg) 3138 if (strg/="") then 3139 strg1=trim(strg) 3140 read(strg1,*)(cell(1,1,mu,2,idisp),mu=1,3) 3141 else 3142 strg1=trim(line) 3143 read(strg1,*)(cell(1,1,mu,2,idisp),mu=1,3) 3144 end if 3145 end if 3146 end if 3147 if ((line(1:19)==char(60)//'/displacement_diff')) then 3148 found2=.false. 3149 end if 3150 end do 3151 end if 3152 end do!end do while displacement 3153 ! In the XML the atom index begin to zero 3154 ! Need to shift for fortran array 3155 atindx(1,1,:,:) = atindx(1,1,:,:) + 1 3156 ! Initialisation of the polynomial_term structure with the values from the 3157 ! previous step 3158 iterm = iterm + 1 3159 call polynomial_term_init(atindx(1,1,:,:),cell(1,1,:,:,:),& 3160 & direction(1,1,:),ndisp,nstrain,terms(1,iterm),& 3161 & power_disp(1,1,:),power_strain(1,1,:),& 3162 & strain(1,1,:),weight(1,1),check=.true.) 3163 end if!end if term 3164 end do!end do while found (coeff) 3165 3166 ! Initialisation of the polynomial_coefficent structure with the values from the 3167 ! previous step 3168 icoeff = icoeff + 1 3169 call polynomial_coeff_init(coefficient(1),nterm,coeffs(icoeff),terms(1,:)) 3170 call polynomial_coeff_getName(name,coeffs(icoeff),symbols,recompute=.true.) 3171 call polynomial_coeff_setName(name,coeffs(icoeff)) 3172 ! Deallocation of the terms array for this coefficient 3173 do jj=1,nterm_max 3174 call polynomial_term_free(terms(1,jj)) 3175 end do 3176 end if!end if line = coefficient 3177 end if!end if ios==0 3178 end do!end do while on file 3179 3180 close(unit=funit) 3181 3182 #endif 3183 ABI_FREE(terms) 3184 ABI_FREE(atindx) 3185 ABI_FREE(coefficient) 3186 ABI_FREE(cell) 3187 ABI_FREE(direction) 3188 ABI_FREE(strain) 3189 ABI_FREE(power_disp) 3190 ABI_FREE(power_strain) 3191 ABI_FREE(weight) 3192 ABI_FREE(symbols) 3193 end if !End if master 3194 3195 !9-MPI BROADCAST 3196 do ii=1,ncoeff 3197 call polynomial_coeff_broadcast(coeffs(ii),master, comm) 3198 end do 3199 3200 !10-checks 3201 3202 !11-debug print 3203 debug = .FALSE. 3204 if(debug)then 3205 do ii=1,ncoeff 3206 do jj=1,coeffs(ii)%nterm 3207 #if defined HAVE_XML 3208 write(200+my_rank,*)"ii,jj,ndisp,nterm",ii,jj,coeffs(ii)%nterm,coeffs(ii)%terms(jj)%ndisp 3209 write(200+my_rank,*)"atindx",coeffs(ii)%terms(jj)%atindx 3210 write(200+my_rank,*)"cell1",coeffs(ii)%terms(jj)%cell(:,1,:) 3211 write(200+my_rank,*)"cell2",coeffs(ii)%terms(jj)%cell(:,2,:) 3212 write(200+my_rank,*)"direction",coeffs(ii)%terms(jj)%direction 3213 write(200+my_rank,*)"power_disp",coeffs(ii)%terms(jj)%power_disp 3214 write(200+my_rank,*)"weight",coeffs(ii)%terms(jj)%weight 3215 #else 3216 write(300+my_rank,*)"ii,jj,ndisp,nterm",ii,jj,coeffs(ii)%nterm,coeffs(ii)%terms(jj)%ndisp 3217 write(300+my_rank,*)"atindx",coeffs(ii)%terms(jj)%atindx 3218 write(300+my_rank,*)"cell1",coeffs(ii)%terms(jj)%cell(:,1,:) 3219 write(300+my_rank,*)"cell2",coeffs(ii)%terms(jj)%cell(:,2,:) 3220 write(300+my_rank,*)"direction",coeffs(ii)%terms(jj)%direction 3221 write(300+my_rank,*)"power_disp",coeffs(ii)%terms(jj)%power_disp 3222 write(300+my_rank,*)"weight",coeffs(ii)%terms(jj)%weight 3223 #endif 3224 end do 3225 end do 3226 #if defined HAVE_XML 3227 close(200+my_rank) 3228 #else 3229 close(300+my_rank) 3230 #endif 3231 end if 3232 3233 !12-Initialisation of eff_pot 3234 call effective_potential_setCoeffs(coeffs,eff_pot,ncoeff) 3235 3236 !13-Deallocation of type 3237 do ii=1,ncoeff 3238 call polynomial_coeff_free(coeffs(ii)) 3239 end do 3240 ABI_FREE(coeffs) 3241 3242 end subroutine coeffs_xml2effpot
m_effective_potential_file/effective_potential_file_getDimCoeff [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_getDimCoeff
FUNCTION
This routine test the xml with polynomial coefficients Return the number of coefficients and the maximum number of displacement/strain
INPUTS
filename = names of the files
OUTPUT
ncoeff = number of coefficient for the polynome nterm(ncoeff) = number terms per coefficient ndisp(nterm,ncoeff) = number displacement per term
SOURCE
670 subroutine effective_potential_file_getDimCoeff(filename,ncoeff,ndisp_max,nterm_max) 671 672 !Arguments ------------------------------------ 673 !scalars 674 character(len=fnlen),intent(in) :: filename 675 integer,intent(out) :: ncoeff,ndisp_max,nterm_max 676 !Local variables------------------------------- 677 !scalar 678 integer :: filetype 679 #ifndef HAVE_XML 680 integer :: count,count2 681 integer :: funit = 1,ios=0 682 logical :: found,found2 683 #endif 684 !arrays 685 #ifndef HAVE_XML 686 character (len=XML_RECL) :: line,readline 687 #endif 688 character(len=500) :: message 689 690 ! ************************************************************************* 691 692 call effective_potential_file_getType(filename,filetype) 693 694 if (filetype==3 .or. filetype==23) then 695 write(message, '(2a)' )' Extraction of the number of coefficient in the XML ',& 696 & trim(filename) 697 call wrtout(std_out,message,'COLL') 698 699 ncoeff = 0 700 nterm_max = 0 701 ndisp_max = 0 702 703 #if defined HAVE_XML 704 ! Read with libxml the number of coefficient 705 call effpot_xml_getDimCoeff(char_f2c(trim(filename)),ncoeff,nterm_max,ndisp_max) 706 #else 707 ! Read by hand 708 ! Start a reading loop 709 found=.false. 710 ncoeff = 0 711 712 if (open_file(filename,message,unit=funit,form="formatted",status="old",& 713 & action="read") /= 0) then 714 ABI_ERROR(message) 715 end if 716 717 ! First parse to know the number of coefficients 718 ios = 0 719 do while (ios == 0) 720 read(funit,'(a)',iostat=ios) readline 721 if(ios == 0)then 722 call rmtabfromline(readline) 723 line=adjustl(readline) 724 ! Need test with char(9) because the old version of XML file 725 ! from old script includes tarbulation at the begining of each line 726 if (line(1:12)==char(60)//'coefficient') then 727 ncoeff=ncoeff+1 728 count = 0 729 found = .false. 730 do while(.not.found) 731 read(funit,'(a)',iostat=ios) readline 732 call rmtabfromline(readline) 733 line=adjustl(readline) 734 if (line(1:5)==char(60)//'term') then 735 count = count +1 736 found2 = .false. 737 count2 = 0 738 do while(.not.found2) 739 read(funit,'(a)',iostat=ios) readline 740 call rmtabfromline(readline) 741 line=adjustl(readline) 742 if (line(1:13)==char(60)//'displacement') then 743 count2 = count2 + 1 744 else if (line(1:7)==char(60)//'strain') then 745 count2 = count2 + 1 746 else if (line(1:6)==char(60)//'/term') then 747 if (count2 > ndisp_max) ndisp_max = count2 748 found2 = .true. 749 else 750 cycle 751 end if 752 end do 753 else if (line(1:13)==char(60)//'/coefficient') then 754 if (count > nterm_max) nterm_max = count 755 found = .true. 756 else 757 cycle 758 end if 759 end do 760 cycle 761 end if 762 end if 763 end do 764 765 close(funit) 766 #endif 767 768 else 769 ! Maybe one day add an other type of file... 770 write(message, '(a,a,a,a)' )& 771 & ' The file ',trim(filename),' is not compatible with multibinit',ch10 772 ABI_ERROR(message) 773 end if 774 775 ! Do some checks 776 if (ncoeff < 1) then 777 write(message, '(5a)' )& 778 & ' Unable to read the number of coeff from ',trim(filename),ch10,& 779 & ' This file is not compatible with multibinit',ch10 780 ABI_ERROR(message) 781 end if 782 783 end subroutine effective_potential_file_getDimCoeff
m_effective_potential_file/effective_potential_file_getDimMD [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_getDimMD
FUNCTION
Read MD FILE (HIST or ASCII) and return the dimensions (natom and nstep)
INPUTS
filename = path of the file
OUTPUT
natom = number of atoms nstep = number of MD steps
SOURCE
899 subroutine effective_potential_file_getDimMD(filename,natom,nstep) 900 901 !Arguments ------------------------------------ 902 !scalars 903 integer,intent(out) :: natom,nstep 904 !arrays 905 character(len=fnlen),intent(in) :: filename 906 !Local variables------------------------------- 907 !scalar 908 integer :: ia,natm_old,natm_new 909 integer :: nenergy,nrprimd 910 integer :: ios=0,ios2=0,ios3=0 911 integer :: unit_md=24 912 logical :: compatible,netcdf 913 #if defined HAVE_NETCDF 914 integer :: natom_id,time_id,xyz_id,six_id 915 integer :: ncid,ncerr 916 character(len=5) :: char_tmp 917 #endif 918 !arrays 919 character (len=10000) :: readline,line 920 character(len=500) :: msg 921 922 ! ************************************************************************* 923 924 natom = 0 925 nstep = 0 926 !try to read netcdf 927 netcdf = .false. 928 #if defined HAVE_NETCDF 929 ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid) 930 if(ncerr == NF90_NOERR) then 931 netcdf = .TRUE. 932 ncerr = nf90_inq_dimid(ncid,"natom",natom_id) 933 if(ncerr /= NF90_NOERR) netcdf = .FALSE. 934 ncerr = nf90_inq_dimid(ncid,"xyz",xyz_id) 935 if(ncerr /= NF90_NOERR) netcdf = .FALSE. 936 ncerr = nf90_inq_dimid(ncid,"time",time_id) 937 if(ncerr /= NF90_NOERR) netcdf = .FALSE. 938 ncerr = nf90_inq_dimid(ncid,"six",six_id) 939 if(ncerr /= NF90_NOERR) netcdf = .FALSE. 940 if(netcdf)then 941 ncerr = nf90_inquire_dimension(ncid,natom_id,char_tmp,natom) 942 NCF_CHECK_MSG(ncerr," inquire dimension ID for natom") 943 ncerr = nf90_inquire_dimension(ncid,time_id,char_tmp,nstep) 944 NCF_CHECK_MSG(ncerr," inquire dimension ID for time") 945 end if 946 end if 947 #endif 948 949 if(.not.netcdf) then 950 ! try to read ASCII file... 951 if (open_file(filename,msg,unit=unit_md,form="formatted",& 952 & status="old",action="read") /= 0) then 953 ABI_ERROR(msg) 954 end if 955 956 ! Start a reading loop in fortran to get the dimension of the file 957 rewind(unit=unit_md) 958 ios = 0 959 nstep = 0 960 nrprimd = 0 961 natm_old= 0 962 natm_new= 0 963 nenergy = 0 964 compatible = .TRUE. 965 966 do while ((ios==0)) 967 ! special treatment of the first step 968 if(nstep==0)then 969 ios2 = 0 970 do while ((ios2==0)) 971 read(unit_md,'(a)',iostat=ios) readline 972 line=adjustl(readline) 973 call elementfromline(line,ia) 974 if (ia==1)then 975 nstep = nstep + 1 976 ios2 = 1 977 end if 978 end do 979 end if 980 read(unit_md,'(a)',iostat=ios) readline 981 if(ios == 0)then 982 line=adjustl(readline) 983 call elementfromline(line,ia) 984 if (ia==1)then 985 nenergy = nenergy + 1 986 nrprimd = 0 987 else if(ia==3)then 988 nrprimd = nrprimd + 1 989 end if 990 if(nrprimd == 3)then 991 ios3 = 0 992 natm_new = 0 993 do while ((ios3==0)) 994 read(unit_md,'(a)',iostat=ios3) readline 995 if(ios3==0)then 996 line=adjustl(readline) 997 call elementfromline(line,ia) 998 if(ia==1)then 999 if(nstep==1) then 1000 natm_old = natm_new 1001 else 1002 if(natm_old /= natm_new) compatible = .FALSE. 1003 end if 1004 ios3 = 1 1005 ios2 = 1 1006 nstep = nstep + 1 1007 end if 1008 if(ia==6)then 1009 natm_new = natm_new + 1 1010 end if 1011 end if!end if ios3 1012 end do 1013 end if ! end if nrprimd 1014 end if! end if os1 1015 end do 1016 1017 natom = natm_new - 1 1018 if(nstep /= nenergy) compatible = .FALSE. 1019 if(natom <= 0) compatible = .FALSE. 1020 if(nstep <= 0) compatible = .FALSE. 1021 1022 if(.not.compatible)then 1023 natom = 0 1024 nstep = 0 1025 end if 1026 end if! end if not netcdf 1027 1028 end subroutine effective_potential_file_getDimMD
m_effective_potential_file/effective_potential_file_getDimStrainCoupling [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_getDimStrainCoupling
FUNCTION
Return the number of nrpt for specific strain coupling from xml system file
INPUTS
filename = names of the files voigt = voigt notation of the strain
OUTPUT
nrpt = number of rpt points
SOURCE
804 subroutine effective_potential_file_getDimStrainCoupling(filename,nrpt,voigt) 805 806 !Arguments ------------------------------------ 807 !scalars 808 character(len=fnlen),intent(in) :: filename 809 integer,intent(in) :: voigt 810 integer,intent(out) :: nrpt 811 !Local variables------------------------------- 812 !scalar 813 #ifndef HAVE_XML 814 integer :: irpt,ivoigt 815 integer :: funit = 1,ios=0 816 logical :: found 817 #endif 818 !arrays 819 #ifndef HAVE_XML 820 character (len=XML_RECL) :: line,readline,strg,strg1 821 character(len=500) :: message 822 #endif 823 824 ! ************************************************************************* 825 826 nrpt = 0 827 828 #if defined HAVE_XML 829 ! Read with libxml the number of coefficient 830 call effpot_xml_getDimStrainCoupling(char_f2c(trim(filename)),nrpt,voigt) 831 #else 832 ! Read by hand 833 ! Start a reading loop 834 found=.false. 835 836 if (open_file(filename,message,unit=funit,form="formatted",status="old",& 837 & action="read") /= 0) then 838 ABI_ERROR(message) 839 end if 840 841 ! First parse to know the number of atoms 842 do while (ios == 0.and.(.not.found)) 843 read(funit,'(a)',iostat=ios) readline 844 if(ios == 0)then 845 call rmtabfromline(readline) 846 line=adjustl(readline) 847 if ((line(1:16)=='<strain_coupling')) then 848 read(funit,'(a)',iostat=ios) readline 849 call rdfromline("voigt",line,strg) 850 strg1=trim(strg) 851 read(strg1,*) ivoigt 852 if (ivoigt == voigt)then 853 irpt = 0 854 do while (.not.found) 855 read(funit,'(a)',iostat=ios) readline 856 call rmtabfromline(readline) 857 line=adjustl(readline) 858 if ((line(1:26)=='<correction_force_constant')) then 859 irpt = irpt + 1 860 cycle 861 end if 862 if ((line(1:17)=='</strain_coupling')) then 863 found = .TRUE. 864 nrpt = irpt 865 cycle 866 end if 867 end do 868 else 869 cycle 870 end if 871 end if 872 end if 873 end do 874 875 close(funit) 876 #endif 877 878 end subroutine effective_potential_file_getDimStrainCoupling
m_effective_potential_file/effective_potential_file_getDimSystem [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_getDimSystem
FUNCTION
This routine test the xml or ddb file Return the number of atoms/ntypat in the unit cell from ddb and xml Return natom/ntypat/nqpt and nrpt if the file is XML file In case of DDB file, you have to run bigbx9 to get nrpt
INPUTS
filename = names of the files comm = MPI communicator
OUTPUT
natom = number of atoms ntypat= number of type of atoms nqpt = number of q points nrpt = number of rpt points
SOURCE
561 subroutine effective_potential_file_getDimSystem(filename,comm,natom,ntypat,nqpt,nrpt) 562 563 !Arguments ------------------------------------ 564 !scalars 565 character(len=fnlen),intent(in) :: filename 566 integer,intent(out) :: natom,ntypat,nqpt,nrpt 567 integer,intent(in) :: comm 568 !arrays 569 570 !Local variables------------------------------- 571 !scalar 572 integer :: filetype 573 ! integer :: dimekb,lmnmax,mband,mtyp,msym,nblok,nkpt,usepaw 574 character(len=500) :: message 575 type(ddb_hdr_type) :: ddb_hdr 576 !arrays 577 ! ************************************************************************* 578 579 natom = 0 580 ntypat= 0 581 nqpt = 0 582 nrpt = 0 583 584 call effective_potential_file_getType(filename,filetype) 585 586 if(filetype==1)then 587 write(message, '(6a)' )ch10,' The file ',trim(filename),ch10,& 588 & ' is DDB file (extraction of the number of atoms)',ch10 589 call wrtout(std_out,message,'COLL') 590 591 write(message, '(8a)' )& 592 & ' The file ',trim(filename),ch10,' is ddb file only the number of atoms is read,',& 593 & 'if you want to predic the number of cell (nrpt)',ch10,' use bigbx9 routines',ch10 594 call wrtout(std_out,message,'COLL') 595 596 call ddb_hdr%open_read(filename,comm,dimonly=1) 597 natom = ddb_hdr%natom 598 ntypat = ddb_hdr%ntypat 599 call ddb_hdr%free() 600 601 ! Must read some value to initialze array (nprt for ifc) 602 ! call bigbx9(inp%brav,dummy_cell,0,1,inp%ngqpt,inp%nqshft,nrpt,ddb%rprim,dummy_rpt) 603 604 else if (filetype==2 .or. filetype==23) then 605 write(message, '(5a)' )ch10,' The file ',trim(filename),& 606 & ' is XML file (extraction of all information)' 607 call wrtout(std_out,message,'COLL') 608 609 call system_getDimFromXML(filename,natom,ntypat,nqpt,nrpt) 610 611 else 612 write(message, '(a,a,a,a)' )& 613 & ' The file ',trim(filename),' is not compatible with multibinit',ch10 614 ABI_ERROR(message) 615 end if 616 617 ! TODO hexu: temporarily disabled. Discuss with alex how to do this properly. 618 ! Do some checks 619 ! if (natom < 1) then 620 ! write(message, '(a,a,a,a,a)' )& 621 !& ' Unable to read the number of atom from ',trim(filename),ch10,& 622 !& 'This file is not compatible with multibinit',ch10 623 ! ABI_ERROR(message) 624 ! end if 625 ! 626 ! if (filetype==2 .or. filetype==23) then 627 ! 628 ! if (natom < 1) then 629 ! write(message, '(a,a,a)' )& 630 !& ' Unable to read the number of atom from ',trim(filename),ch10 631 ! ABI_ERROR(message) 632 ! end if 633 ! 634 ! if (nrpt < 1) then 635 ! write(message, '(a,a,a)' )& 636 !& ' Unable to read the number of rpt points ',trim(filename),ch10 637 ! ABI_ERROR(message) 638 ! end if 639 ! 640 ! if (ntypat < 1) then 641 ! write(message, '(a,a,a)' )& 642 !& ' Unable to read the number of type of atoms ',trim(filename),ch10 643 ! ABI_ERROR(message) 644 ! end if 645 ! 646 ! end if 647 648 end subroutine effective_potential_file_getDimSystem
m_effective_potential_file/effective_potential_file_getType [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_getType
FUNCTION
This routine test the xml or ddb file
INPUTS
filename = names of the files
OUTPUT
type_file = 0 no type found
1 DDB file
2 XML file the system definition and harmonic part
3 XML file with polynomial coefficients
23 XML file with both system definition and polynomial coefficients
40 NetCDF file with history of MD or snapshot
41 ASCII file with history of MD or snapshot
SOURCE
427 subroutine effective_potential_file_getType(filename,filetype) 428 429 !Arguments ------------------------------------ 430 !scalars 431 character(len=fnlen),intent(in) :: filename 432 integer, intent(out) :: filetype 433 !arrays 434 !Local variables------------------------------- 435 !scalar 436 integer :: natom,nstep 437 integer :: ddbun = 666,ios=0 438 character(len=500) :: message 439 character (len=1000) :: line,readline 440 #if defined HAVE_NETCDF 441 integer :: natom_id,time_id,xyz_id,six_id,ddb_version 442 integer :: ncid,ncerr 443 logical :: md_file 444 #endif 445 446 !arrays 447 ! ************************************************************************* 448 449 filetype = 0 450 451 ddbun = get_unit() 452 453 if (open_file(filename,message,unit=ddbun,form="formatted",status="old",action="read") /= 0) then 454 ABI_ERROR(message) 455 end if 456 457 !Check if the file is a XML file or a DDB and in this case, store the DDB code. 458 ios = 0 459 do while ((ios==0)) 460 read(ddbun,'(a)',iostat=ios) readline 461 call rmtabfromline(readline) 462 line=adjustl(readline) 463 if(line(3:13)=="xml version") then 464 do while ((ios==0)) 465 read(ddbun,'(a)',iostat=ios) readline 466 call rmtabfromline(readline) 467 line=adjustl(readline) 468 if(line(1:16)==char(60)//"Heff_definition".or.& 469 & line(1:17)==char(60)//"Terms_definition")then 470 filetype = 3 471 ios = -1 472 end if 473 if(line(1:18)==char(60)//"System_definition") then 474 filetype = 2 475 do while ((ios==0)) 476 read(ddbun,'(a)',iostat=ios) readline 477 call rmtabfromline(readline) 478 line=adjustl(readline) 479 if(line(1:16)==char(60)//"Heff_definition".or.& 480 & line(1:17)==char(60)//"Terms_definition")then 481 filetype = 23 482 ios = -1 483 end if 484 end do 485 end if 486 end do 487 else if(line(6:24)=="DERIVATIVE DATABASE") then 488 filetype = 1 489 ios = -1 490 end if 491 end do 492 close(ddbun) 493 494 if(filetype/=0) return 495 496 !try to read netcdf HIST file 497 #if defined HAVE_NETCDF 498 ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid) 499 if(ncerr == NF90_NOERR) then 500 md_file = .TRUE. 501 ncerr = nf90_inq_dimid(ncid,"natom",natom_id) 502 if(ncerr /= NF90_NOERR) md_file = .FALSE. 503 ncerr = nf90_inq_dimid(ncid,"xyz",xyz_id) 504 if(ncerr /= NF90_NOERR) md_file = .FALSE. 505 ncerr = nf90_inq_dimid(ncid,"time",time_id) 506 if(ncerr /= NF90_NOERR) md_file = .FALSE. 507 ncerr = nf90_inq_dimid(ncid,"six",six_id) 508 if(ncerr /= NF90_NOERR) md_file = .FALSE. 509 if (md_file) then 510 filetype = 40 511 return 512 end if 513 end if 514 ncerr = nf90_close(ncid) 515 #endif 516 517 if(filetype/=0) return 518 519 ! Try to read netcdf DDB file 520 #if defined HAVE_NETCDF 521 ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid) 522 if(ncerr==NF90_NOERR) then 523 ncerr = nf90_get_var(ncid, nctk_idname(ncid, 'ddb_version'), ddb_version) 524 if (ncerr==NF90_NOERR) then 525 filetype = 1 526 return 527 end if 528 end if 529 #endif 530 531 !Try to get the dim of MD ASCII file 532 call effective_potential_file_getDimMD(filename,natom,nstep) 533 if(natom /= 0 .and. nstep/=0) filetype = 41 534 535 end subroutine effective_potential_file_getType
m_effective_potential_file/effective_potential_file_mapHistToRef [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_mapHistToRef
FUNCTION
Generate the supercell in the effective potential according to the size of the supercell in the hist file Check if the hist file match to reference supercell in the effective potential If not, the hist file is reordering
INPUTS
eff_pot<type(effective_potential)> = effective potential hist<type(abihist)> = The history of the MD comm = MPI communicator
OUTPUT
hist<type(abihist)> = The history of the MD
SOURCE
3372 subroutine effective_potential_file_mapHistToRef(eff_pot,hist,comm,iatfix,verbose,sc_size) 3373 3374 !Arguments ------------------------------------ 3375 !scalars 3376 integer,intent(in) :: comm 3377 logical,optional,intent(in) :: verbose 3378 !arrays 3379 type(effective_potential_type),intent(inout) :: eff_pot 3380 type(abihist),intent(inout) :: hist 3381 integer,optional,allocatable,intent(inout) :: iatfix(:,:) 3382 integer,optional,intent(in) :: sc_size(3) 3383 !Local variables------------------------------- 3384 !scalar 3385 integer :: factE_hist,ia,ib,ii,jj,natom_hist,ncells,nstep_hist 3386 real(dp):: factor,ratio 3387 logical :: revelant_factor,need_map,need_verbose,need_fixmap 3388 !arrays 3389 real(dp) :: rprimd_hist(3,3),rprimd_ref(3,3) 3390 integer :: ncell(3),scale_cell(3) 3391 integer,allocatable :: shift(:,:),iatfix_tmp(:,:) 3392 integer,allocatable :: list_map(:) !blkval(:), 3393 real(dp),allocatable :: xred_ref(:,:) ! xred_hist(:,:), 3394 real(dp),allocatable :: list_dist(:),list_reddist(:,:),list_absdist(:,:) 3395 character(len=500) :: msg 3396 type(abihist) :: hist_tmp 3397 ! ************************************************************************* 3398 3399 !Set optional values 3400 need_verbose = .false. 3401 need_fixmap = .FALSE. 3402 if (present(verbose)) need_verbose = verbose 3403 if (present(iatfix)) need_fixmap = .TRUE. 3404 3405 natom_hist = size(hist%xred,2) 3406 nstep_hist = size(hist%xred,3) 3407 3408 ! Try to set the supercell according to the hist file 3409 rprimd_ref(:,:) = eff_pot%crystal%rprimd 3410 rprimd_hist(:,:) = hist%rprimd(:,:,1) 3411 3412 if(present(sc_size))then 3413 ncell(:) = sc_size 3414 else 3415 do ia=1,3 3416 scale_cell(:) = 0 3417 do ii=1,3 3418 if(abs(rprimd_ref(ii,ia)) > tol10)then 3419 scale_cell(ii) = nint(rprimd_hist(ii,ia) / rprimd_ref(ii,ia)) 3420 end if 3421 end do 3422 ! Check if the factor for the supercell is revelant 3423 revelant_factor = .TRUE. 3424 do ii=1,3 3425 if(abs(scale_cell(ii)) < tol10) cycle 3426 factor = abs(scale_cell(ii)) 3427 do jj=ii,3 3428 if(abs(scale_cell(jj)) < tol10) cycle 3429 if(abs(abs(scale_cell(ii))-abs(scale_cell(jj))) > tol10) revelant_factor = .FALSE. 3430 end do 3431 end do 3432 if(.not.revelant_factor)then 3433 write(msg, '(3a)' )& 3434 & 'unable to map the hist file ',ch10,& 3435 & 'Action: check/change your MD file' 3436 ABI_ERROR(msg) 3437 else 3438 ncell(ia) = int(factor) 3439 end if 3440 end do 3441 end if 3442 3443 ncells = product(ncell) 3444 3445 !Check if the energy stored in the hist is revelant, sometimes some MD files gives 3446 !the energy of the unit cell... This is not suppose to happen... But just in case... 3447 do ii=1,nstep_hist 3448 if(abs(eff_pot%energy)>tol12)then 3449 ratio=hist%etot(ii) / eff_pot%energy 3450 if(abs(ratio)<real(huge(factE_hist))*half)then 3451 factE_hist = nint(ratio) 3452 if(factE_hist == 1) then 3453 ! In this case we mutiply the energy of the hist by the number of cell 3454 hist%etot(ii) = hist%etot(ii) * ncells 3455 end if 3456 if(factE_hist /=1 .and. factE_hist /= ncells)then 3457 write(msg, '(4a,I0,a,I0,2a,I0,3a,I0,3a)' )ch10,& 3458 & ' --- !WARNING',ch10,& 3459 & ' The energy of the history step ',ii,' seems to be with multiplicity of ',factE_hist,ch10,& 3460 & ' However, the multiplicity of the cell is ',ncells,'.',ch10,& 3461 & ' Please check the energy of the step ',ii,ch10,& 3462 & ' ---',ch10 3463 if(need_verbose) call wrtout(std_out,msg,'COLL') 3464 endif 3465 else 3466 write(msg, '(4a,i0,3a,es16.6,5a)' )ch10,& 3467 & ' --- !WARNING',ch10,& 3468 & ' The energy of the history step ',ii,' is apparently not initialized.',ch10,& 3469 & ' Its current value is',hist%etot(ii),ch10,& 3470 & ' This does not allow to perform checking on the multiplicity of the cell ',ch10,& 3471 & ' ---',ch10 3472 if(need_verbose) call wrtout(std_out,msg,'COLL') 3473 end if 3474 end if 3475 end do 3476 3477 3478 !Set the new supercell datatype into the effective potential reference 3479 call effective_potential_setSupercell(eff_pot,comm,ncell) 3480 3481 !allocation 3482 ABI_MALLOC(shift,(3,natom_hist)) 3483 ABI_MALLOC(list_map,(natom_hist)) 3484 ABI_MALLOC(list_reddist,(3,natom_hist)) 3485 ABI_MALLOC(list_absdist,(3,natom_hist)) 3486 ABI_MALLOC(list_dist,(natom_hist)) 3487 ABI_MALLOC(xred_ref,(3,natom_hist)) 3488 3489 !Putting maping list to zero 3490 list_map = 0 3491 3492 !Fill xcart_ref/hist and xred_ref/hist 3493 3494 call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,& 3495 & eff_pot%supercell%xcart,xred_ref) ! Get xred_ref 3496 3497 3498 do ia=1,natom_hist !Loop over all reference atoms 3499 ! Put temporary lists to zero 3500 list_reddist = 0 3501 list_absdist = 0 3502 list_dist = 0 3503 shift = 0 3504 do ib=1,natom_hist !Loop over all atoms of distorted structure 3505 !Calculate list of reduced distance between reference atom ia and all others 3506 list_reddist(:,ib) = hist%xred(:,ib,1) - xred_ref(:,ia) 3507 !If the distorted atom is further away than half the unit cell shift it. 3508 if(list_reddist(1,ib) > 0.5)then 3509 list_reddist(1,ib) = 1 - list_reddist(1,ib) 3510 shift(1,ib) = -1 3511 end if 3512 if(list_reddist(2,ib) > 0.5)then 3513 list_reddist(2,ib) = 1 - list_reddist(2,ib) 3514 shift(2,ib) = -1 3515 end if 3516 if(list_reddist(3,ib) > 0.5)then 3517 list_reddist(3,ib) = 1 - list_reddist(3,ib) 3518 shift(3,ib) = -1 3519 end if 3520 if(list_reddist(1,ib) < -0.5)then 3521 list_reddist(1,ib) = -1 - list_reddist(1,ib) 3522 shift(1,ib) = 1 3523 end if 3524 if(list_reddist(2,ib) < -0.5)then 3525 list_reddist(2,ib) = -1 - list_reddist(2,ib) 3526 shift(2,ib) = 1 3527 end if 3528 if(list_reddist(3,ib) < -0.5)then 3529 list_reddist(3,ib) = -1 - list_reddist(3,ib) 3530 shift(3,ib) = 1 3531 end if 3532 list_absdist(1,ib) = (rprimd_hist(1,1)+rprimd_hist(2,1)+rprimd_hist(3,1))*list_reddist(1,ib) 3533 list_absdist(2,ib) = (rprimd_hist(1,2)+rprimd_hist(2,2)+rprimd_hist(3,2))*list_reddist(2,ib) 3534 list_absdist(3,ib) = (rprimd_hist(1,3)+rprimd_hist(2,3)+rprimd_hist(3,3))*list_reddist(3,ib) 3535 list_dist(ib) = sqrt(abs(list_absdist(1,ib))**2 + abs(list_absdist(2,ib))**2 + abs(list_absdist(3,ib))**2 ) 3536 end do !ib 3537 !find the closest atom ib 3538 list_map(ia) = minloc(list_dist,DIM=1) 3539 !If the closest atom ib was shifted, apply and store the shift 3540 if(any(shift(:,list_map(ia)) /= 0))then 3541 hist%xred(1,list_map(ia),:)= hist%xred(1,list_map(ia),:) + 1*shift(1,list_map(ia)) 3542 hist%xred(2,list_map(ia),:)= hist%xred(2,list_map(ia),:) + 1*shift(2,list_map(ia)) 3543 hist%xred(3,list_map(ia),:)= hist%xred(3,list_map(ia),:) + 1*shift(3,list_map(ia)) 3544 end if 3545 !TEST MS 3546 !write(*,*) 'Atom', ia,' of reference is matche with', list_map(ia) 3547 !write(*,*) 'xred_ref(',xred_ref(:,ia),'), xred_hist(',hist%(:,list_map(ia)),')' 3548 end do ! ia 3549 3550 if(need_verbose) then 3551 write(msg,'(2a,I3,a,I3,a,I3)') ch10,& 3552 & ' The size of the supercell for the fit is ',ncell(1),' ',ncell(2),' ',ncell(3) 3553 call wrtout(std_out,msg,'COLL') 3554 call wrtout(ab_out,msg,'COLL') 3555 end if 3556 3557 3558 if(any(list_map(:)==0))then 3559 write(msg, '(5a)' )& 3560 & 'Unable to map the molecular dynamic file ',ch10,& 3561 & 'on the reference supercell structure',ch10,& 3562 & 'Action: change the MD file' 3563 ABI_ERROR(msg) 3564 end if 3565 3566 need_map = .FALSE. 3567 do ia=1,natom_hist 3568 if(list_map(ia) /= ia) need_map = .TRUE. 3569 end do 3570 if(need_map)then 3571 if(need_verbose) then 3572 write(msg, '(11a)' )ch10,& 3573 & ' --- !WARNING',ch10,& 3574 & ' The ordering of the atoms in the _HIST.nc file is different,',ch10,& 3575 & ' of the one built by multibinit. The _HIST.nc file will be mapped,',ch10,& 3576 & ' to the ordering of multibinit.',ch10,& 3577 & ' ---',ch10 3578 call wrtout(ab_out,msg,'COLL') 3579 call wrtout(std_out,msg,'COLL') 3580 end if 3581 3582 ! Allocate hist datatype 3583 call abihist_init(hist_tmp,natom_hist,nstep_hist,.false.,.false.) 3584 ! copy all the information 3585 do ia=1,nstep_hist 3586 hist%ihist = ia 3587 hist_tmp%ihist = ia 3588 call abihist_copy(hist,hist_tmp) 3589 end do 3590 hist_tmp%mxhist = nstep_hist 3591 3592 ! reoder array 3593 do ia=1,natom_hist 3594 hist_tmp%xred(:,ia,:) = hist%xred(: ,list_map(ia),:) 3595 hist_tmp%fcart(:,ia,:) = hist%fcart(:,list_map(ia),:) 3596 hist_tmp%vel(:,ia,:) = hist%vel(:,list_map(ia),:) 3597 end do 3598 3599 ! free the old hist and reinit 3600 call abihist_free(hist) 3601 call abihist_init(hist,natom_hist,nstep_hist,.false.,.false.) 3602 ! copy the temporary hist into output 3603 do ia=1,nstep_hist 3604 hist%ihist = ia 3605 hist_tmp%ihist = ia 3606 call abihist_copy(hist_tmp,hist) 3607 end do 3608 hist_tmp%mxhist = nstep_hist 3609 call abihist_free(hist_tmp) 3610 3611 !map also fixes if present 3612 if(need_fixmap)then 3613 ABI_MALLOC(iatfix_tmp,(3,natom_hist)) 3614 do ia=1,natom_hist 3615 iatfix_tmp(:,ia) = iatfix(:,list_map(ia)) 3616 end do 3617 iatfix = iatfix_tmp 3618 ABI_FREE(iatfix_tmp) 3619 end if 3620 end if !need map 3621 3622 !deallocation 3623 ABI_FREE(shift) 3624 ABI_FREE(list_map) 3625 ABI_FREE(list_dist) 3626 ABI_FREE(list_reddist) 3627 ABI_FREE(list_absdist) 3628 ABI_FREE(xred_ref) 3629 end subroutine effective_potential_file_mapHistToRef
m_effective_potential_file/effective_potential_file_readMDfile [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
effective_potential_file_readMDfile
FUNCTION
Read MD FILE (HIST or ASCII)
INPUTS
filename = path of the file option,optional = 0 (default), the stress is printed in the MD File 1, the force on the cell is printed in the MD File (-1 * stress), in this case, we multiply the stress by -1 in order to get the stresse
OUTPUT
hist<type(abihist)> = datatype with the history of the MD
SOURCE
3264 subroutine effective_potential_file_readMDfile(filename,hist,option) 3265 3266 !Arguments ------------------------------------ 3267 !scalars 3268 integer,optional :: option 3269 !arrays 3270 type(abihist),intent(inout) :: hist 3271 character(len=fnlen),intent(in) :: filename 3272 !Local variables------------------------------- 3273 !scalar 3274 integer :: ia,ii,mu,nu,natom,nstep,type,option_in 3275 integer :: ios=0, unit_md=24 3276 !arrays 3277 character (len=10000) :: readline,line 3278 real(dp) :: tmp(6) 3279 real(dp),allocatable :: xcart(:,:) 3280 3281 ! ************************************************************************* 3282 3283 call effective_potential_file_getType(filename,type) 3284 3285 option_in = 0 3286 if(present(option))then 3287 option_in = option 3288 end if 3289 if(type==40)then 3290 ! Netcdf type 3291 call read_md_hist(filename,hist,.FALSE.,.FALSE.,.FALSE.) 3292 3293 else if(type==41)then 3294 3295 ! ASCII file 3296 call effective_potential_file_getDimMD(filename,natom,nstep) 3297 3298 ii = 1 3299 ios = 0 3300 3301 ABI_MALLOC(xcart,(3,natom)) 3302 call abihist_free(hist) 3303 call abihist_init(hist,natom,nstep,.FALSE.,.FALSE.) 3304 3305 ! Start a reading loop in fortran 3306 rewind(unit=unit_md) 3307 do while ((ios==0).and.ii<=nstep) 3308 read(unit_md,'(a)',iostat=ios) readline 3309 read(unit_md,'(a)',iostat=ios) readline 3310 line=adjustl(readline) 3311 read(line,*) hist%etot(ii) 3312 hist%etot(ii) = hist%etot(ii) 3313 do mu=1,3 3314 read(unit_md,'(a)',iostat=ios) readline 3315 line=adjustl(readline) 3316 read(line,*) (hist%rprimd(nu,mu,ii),nu=1,3) 3317 end do 3318 do ia=1,natom 3319 read(unit_md,'(a)',iostat=ios) readline 3320 line=adjustl(readline) 3321 read(line,*) (tmp(mu),mu=1,6) 3322 xcart(:,ia) = tmp(1:3) 3323 hist%fcart(:,ia,ii) = tmp(4:6) 3324 end do 3325 call xcart2xred(natom,hist%rprimd(:,:,ii),xcart(:,:),hist%xred(:,:,ii)) 3326 read(unit_md,'(a)',iostat=ios) readline 3327 line=adjustl(readline) 3328 read(line,*) (hist%strten(mu,ii),mu=1,6) 3329 ii = ii + 1 3330 end do 3331 do ii=1,nstep 3332 do mu=1,3 3333 hist%acell(mu,:) = hist%rprimd(mu,mu,ii) 3334 end do 3335 end do 3336 close(unit_md) 3337 ABI_FREE(xcart) 3338 3339 end if!end if type 3340 3341 if((type==40 .or. type==41).and.option == 1)then 3342 ! multiply by -1 if the current strten -1*stress, we need only stress... 3343 hist%strten(:,:) = -1 * hist%strten(:,:) 3344 end if 3345 3346 3347 3348 end subroutine effective_potential_file_readMDfile
m_effective_potential_file/elementfromline [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
elementfromline
FUNCTION
Read the number of element of a line
INPUTS
line= string from which the data are read
OUTPUT
nelement = number of element in the line
SOURCE
3708 subroutine elementfromline(line,nelement) 3709 3710 !Arguments --------------------------------------------- 3711 character(len=*), intent(in) :: line 3712 integer, intent(out) :: nelement 3713 !Local variables --------------------------------------- 3714 integer :: ii,n 3715 logical :: element 3716 ! ********************************************************************* 3717 3718 !Set the output 3719 nelement = 0 3720 n = len_trim(line) 3721 element = .false. 3722 do ii=1,n 3723 if(.not.element.and.line(ii:ii) /="") then 3724 element=.true. 3725 else 3726 if((element.and.line(ii:ii) =="")) then 3727 element=.false. 3728 nelement = nelement + 1 3729 end if 3730 end if 3731 if((element.and.ii==n)) nelement = nelement + 1 3732 end do 3733 3734 end subroutine elementfromline
m_effective_potential_file/rdfromline [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
rdfromline
FUNCTION
Read the value of a keyword from a XML line Same function than m_pawxmlps/paw_rdfromline.F90
INPUTS
keyword= keyword which value has to be read line= string from which the data are read (line from a XML)
OUTPUT
output= (string) value of the keyword
SOURCE
3754 subroutine rdfromline(keyword,line,output) 3755 3756 !Arguments --------------------------------------------- 3757 character(len=*), intent(in) :: keyword,line 3758 character(len=*), intent(out) :: output 3759 !Local variables --------------------------------------- 3760 character(len=len(line)) :: temp 3761 integer :: pos,pos2 3762 3763 ! ********************************************************************* 3764 3765 output="" 3766 pos=index(line,trim(keyword)) 3767 if (pos>0) then 3768 temp=line(pos+len_trim(keyword):len_trim(line)) 3769 pos=index(temp,char(34)) 3770 if (pos>0) then 3771 pos2=index(temp(pos+1:len_trim(temp)),char(34)) 3772 if (pos2>0) then 3773 output=temp(pos+1:pos+pos2-1) 3774 end if 3775 end if 3776 end if 3777 3778 end subroutine rdfromline
m_effective_potential_file/rdfromline_value [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
rdfromline
FUNCTION
Read the value of a keyword from a XML line
INPUTS
keyword= keyword which value has to be read line= string from which the data are read (line from a XML)
OUTPUT
output= (string) value of the keyword
SOURCE
3833 subroutine rdfromline_value(keyword,line,output) 3834 3835 !Arguments --------------------------------------------- 3836 character(len=*), intent(in) :: keyword,line 3837 character(len=*), intent(out) :: output 3838 !Local variables --------------------------------------- 3839 character(len=len(line)) :: temp 3840 integer :: pos,pos2 3841 3842 ! ********************************************************************* 3843 3844 output="" 3845 pos=index(line,trim(keyword)) 3846 if (pos==2) then 3847 pos=pos+len_trim(keyword) 3848 pos=pos+index(line(pos:len_trim(line)),char(62)) 3849 temp=line(pos:len_trim(line)) 3850 pos2=index(temp,char(60)) 3851 if (pos2>0) then 3852 output=line(pos:pos+pos2-2) 3853 else 3854 output=line(pos:len_trim(line)) 3855 end if 3856 else 3857 if(pos>2)then 3858 output=line(1:pos-3) 3859 end if 3860 end if 3861 end subroutine rdfromline_value
m_effective_potential_file/rmtabfromline [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
rmtabfromline
FUNCTION
Read remove tab from the begining of line
INPUTS
line= string from which the data are read (line from a XML)
OUTPUT
output= line without tab
SOURCE
3797 recursive subroutine rmtabfromline(line) 3798 3799 !Arguments --------------------------------------------- 3800 character(len=*), intent(inout) :: line 3801 !Local variables --------------------------------------- 3802 integer :: pos 3803 3804 ! ********************************************************************* 3805 3806 pos=index(line,char(9)) 3807 if (pos==1) then 3808 line = line(2:len_trim(line))//" " 3809 call rmtabfromline(line) 3810 end if 3811 3812 end subroutine rmtabfromline
m_effective_potential_file/system_ddb2effpot [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
system_ddb2effpot
FUNCTION
Transfert ddb into effective potential structure. Also calculate the IFC
INPUTS
crytal<type(crystal_t)> = datatype with all the information for the crystal ddb<type(ddb_type)> = datatype with the ddb inp<type(multibinit_dtset_type)> = datatype with the input variables of multibinit comm = MPI communicator
OUTPUT
effective_potantial<type(effective_potential_type)> = effective_potential datatype to be initialized
SOURCE
2162 subroutine system_ddb2effpot(crystal,ddb, effective_potential,inp,comm) 2163 2164 use m_dynmat 2165 2166 use m_copy, only : alloc_copy 2167 use m_crystal, only : crystal_t 2168 use m_multibinit_dataset, only : multibinit_dtset_type 2169 2170 !Arguments ------------------------------------ 2171 !scalars 2172 integer,intent(in) :: comm 2173 !arrays 2174 type(ddb_type),intent(inout) :: ddb 2175 type(effective_potential_type), intent(inout) :: effective_potential 2176 type(crystal_t),intent(in) :: crystal 2177 type(multibinit_dtset_type),intent(in) :: inp 2178 2179 !Local variables------------------------------- 2180 !scalar 2181 real(dp):: wcount1,wcount2 2182 integer :: chneut,i1,i2,i3,ia,ib,iblok,idir1,idir2,ierr,ii,ipert1,iphl1 2183 integer :: ipert2,irpt,irpt2,ivarA,ivarB,max1,max2,max3,min1,min2,min3 2184 integer :: msize,mpert,natom,nblok,nrpt_new,nrpt_new2,rftyp,selectz 2185 integer :: my_rank,nproc,prt_internalstr 2186 logical :: iam_master 2187 integer,parameter :: master=0 2188 integer :: nptsym,nsym 2189 integer :: msym = 384, use_inversion = 1, space_group 2190 real(dp):: max_phfq,tolsym = tol8 2191 !arrays 2192 integer :: bravais(11),cell_number(3),cell2(3) 2193 integer :: shift(3),rfelfd(4),rfphon(4),rfstrs(4) 2194 integer,allocatable :: cell_red(:,:) 2195 real(dp):: dielt(3,3),elast_clamped(6,6),fact 2196 real(dp):: red(3,3),qphnrm(3),qphon(3,3) 2197 real(dp),allocatable :: blkval(:,:,:,:,:,:),d2asr(:,:,:,:,:) 2198 real(dp),allocatable :: instrain(:,:),zeff(:,:,:),qdrp_cart(:,:,:,:) 2199 real(dp),pointer :: atmfrc_red(:,:,:,:,:),wghatm_red(:,:,:) 2200 character(len=500) :: message 2201 type(asrq0_t) :: asrq0 2202 type(ifc_type) :: ifc 2203 real(dp),allocatable :: d2cart(:,:,:,:,:),displ(:) 2204 real(dp),allocatable :: eigval(:,:),eigvec(:,:,:,:,:),phfrq(:) 2205 real(dp),allocatable :: spinat(:,:),tnons(:,:) 2206 integer,allocatable :: symrel(:,:,:),symafm(:),ptsymrel(:,:,:) 2207 2208 ! ************************************************************************* 2209 2210 !0 MPI variables 2211 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2212 iam_master=.FALSE. 2213 iam_master = (my_rank == master) 2214 2215 !Free the eff_pot before filling 2216 call effective_potential_free(effective_potential) 2217 2218 !Initialisation of usefull values 2219 natom = ddb%natom 2220 nblok = ddb%nblok 2221 mpert= ddb%mpert 2222 msize=3*mpert*3*mpert; 2223 2224 !Tranfert the ddb into usable array (ipert and idir format like in abinit) 2225 ABI_MALLOC(blkval,(2,3,mpert,3,mpert,nblok)) 2226 blkval = 0 2227 blkval = reshape(ddb%val,(/2,3,mpert,3,mpert,nblok/)) 2228 2229 !********************************************************************** 2230 ! Transfert crystal values 2231 !********************************************************************** 2232 ! Re-generate symmetry operations from the lattice and atomic coordinates 2233 ABI_MALLOC(spinat,(3,natom)) 2234 ABI_MALLOC(ptsymrel,(3,3,msym)) 2235 ABI_MALLOC(symafm,(msym)) 2236 ABI_MALLOC(symrel,(3,3,msym)) 2237 ABI_MALLOC(tnons,(3,msym)) 2238 spinat = zero; symrel = 0; symafm = 0; tnons = zero ; space_group = 0; 2239 call symlatt(bravais,std_out,msym,nptsym,ptsymrel,crystal%rprimd,tolsym) 2240 call symfind(crystal%gprimd,msym,crystal%natom,nptsym,0,nsym,& 2241 & 0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,& 2242 & crystal%typat,use_inversion,crystal%xred) 2243 if(crystal%nsym/=nsym)then 2244 write(message,'(4a,I0,3a,I0,3a)') ch10,& 2245 & ' --- !WARNING:',ch10,& 2246 & ' There is ',nsym,' found for the crystal',ch10,& 2247 & ' but ',crystal%nsym,' found in the DDB',ch10,& 2248 & ' ---' 2249 call wrtout(std_out,message,'COLL') 2250 end if 2251 call crystal_init(ddb%amu,effective_potential%crystal,& 2252 & space_group,crystal%natom,crystal%npsp,& 2253 & crystal%ntypat,nsym,crystal%rprimd,& 2254 & crystal%typat,crystal%xred,crystal%zion,& 2255 & crystal%znucl,crystal%timrev,crystal%use_antiferro,& 2256 & .FALSE.,crystal%title,& 2257 & symrel=symrel,tnons=tnons,& 2258 & symafm=symafm) 2259 2260 ABI_FREE(spinat) 2261 ABI_FREE(ptsymrel) 2262 ABI_FREE(symafm) 2263 ABI_FREE(symrel) 2264 ABI_FREE(tnons) 2265 2266 !********************************************************************** 2267 ! Transfert energy from input file 2268 !********************************************************************** 2269 write(message, '(2a,(80a),6a)') ch10,('=',ii=1,80),ch10,ch10,& 2270 & ' Extraction of the energy of the structure (unit: Hartree)',ch10 2271 call wrtout(std_out,message,'COLL') 2272 call wrtout(ab_out,message,'COLL') 2273 if (ddb%get_etotal(effective_potential%energy) == 0) then 2274 if(abs(inp%energy_reference) < tol16)then 2275 write(message,'(5a)')& 2276 & ' Warning : Energy of the reference structure is not specify in',& 2277 & ' the input file.',ch10,' Energy will set to zero',ch10 2278 call wrtout(std_out,message,'COLL') 2279 effective_potential%energy = zero 2280 else 2281 effective_potential%energy = inp%energy_reference 2282 end if 2283 else 2284 if(abs(inp%energy_reference) > tol16)then 2285 write(message,'(6a)')& 2286 & ' Warning : Energy of the reference structure is specify in',& 2287 & ' the input file.',ch10,' and in the DDB.',& 2288 & ' The value of the energy is set with the value from the input file',ch10 2289 call wrtout(std_out,message,'COLL') 2290 effective_potential%energy = inp%energy_reference 2291 end if 2292 end if 2293 write(message,'(a,es25.12)') ' Energy = ',& 2294 & effective_potential%energy 2295 call wrtout(std_out,message,'COLL') 2296 call wrtout(ab_out,message,'COLL') 2297 2298 !********************************************************************** 2299 ! Dielectric Tensor and Effective Charges 2300 !********************************************************************** 2301 ABI_MALLOC(zeff,(3,3,natom)) 2302 ABI_MALLOC(qdrp_cart,(3,3,3,natom)) 2303 ABI_MALLOC(effective_potential%harmonics_terms%zeff,(3,3,natom)) 2304 2305 rftyp = 1 ! Blocks obtained by a non-stationary formulation. 2306 chneut = 1 ! The ASR for effective charges is imposed 2307 selectz = 0 ! No selection of some parts of the effective charge tensor 2308 iblok = ddb%get_dielt_zeff(crystal,rftyp,chneut,selectz,dielt,zeff) 2309 qdrp_cart = zero 2310 if (iblok /=0 .and. maxval(abs(dielt)) < 10000) then 2311 effective_potential%harmonics_terms%epsilon_inf = dielt 2312 effective_potential%harmonics_terms%zeff = zeff 2313 else 2314 effective_potential%harmonics_terms%epsilon_inf(1,1) = one 2315 effective_potential%harmonics_terms%epsilon_inf(2,2) = one 2316 effective_potential%harmonics_terms%epsilon_inf(3,3) = one 2317 effective_potential%harmonics_terms%zeff = zero 2318 end if 2319 2320 !********************************************************************** 2321 ! Look after the blok no. that contains the stress tensor 2322 !********************************************************************** 2323 write(message, '(a,a,(80a),a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,& 2324 & ' Extraction of the stress tensor (unit: GPa) and forces (unit: Ha/bohr)' 2325 call wrtout(std_out,message,'COLL') 2326 call wrtout(ab_out,message,'COLL') 2327 2328 ABI_MALLOC(effective_potential%fcart,(3,natom)) 2329 effective_potential%fcart = zero 2330 effective_potential%strten = zero 2331 2332 qphon(:,1)=zero 2333 qphnrm(1)=zero 2334 rfphon(1:2)=0 2335 rfelfd(1:2)=0 2336 rfstrs(1:2)=0 2337 rftyp=4 2338 2339 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 2340 2341 if (iblok /=0) then 2342 if(any(abs(inp%strten_reference)>tol16))then 2343 write(message,'(10a)') ch10,& 2344 & ' --- !WARNING:',ch10,& 2345 & ' The stress tensor of the reference structure is specify in the',ch10,& 2346 & ' input file and in the DDB. The value of the stress tensor is set',ch10,& 2347 & ' with the value from the input file',ch10,& 2348 & ' ---' 2349 call wrtout(std_out,message,'COLL') 2350 call wrtout(ab_out,message,'COLL') 2351 effective_potential%strten = inp%strten_reference 2352 else 2353 ! firts give the corect stress values store in hartree 2354 ! diagonal parts 2355 effective_potential%strten(1)=blkval(1,1,natom+3,1,1,iblok) * crystal%ucvol 2356 effective_potential%strten(2)=blkval(1,2,natom+3,1,1,iblok) * crystal%ucvol 2357 effective_potential%strten(3)=blkval(1,3,natom+3,1,1,iblok) * crystal%ucvol 2358 ! the shear parts 2359 effective_potential%strten(4)=blkval(1,1,natom+4,1,1,iblok) * crystal%ucvol 2360 effective_potential%strten(5)=blkval(1,2,natom+4,1,1,iblok) * crystal%ucvol 2361 effective_potential%strten(6)=blkval(1,3,natom+4,1,1,iblok) * crystal%ucvol 2362 end if 2363 ! Get forces 2364 effective_potential%fcart(:,1:natom) = blkval(1,:,1:natom,1,1,iblok) 2365 else 2366 if(all(abs(inp%strten_reference(:))<tol16))then 2367 write(message,'(8a)') ch10,& 2368 & ' --- !WARNING:',ch10,& 2369 & ' The stress tensor of the reference structure is not specify',ch10,& 2370 & ' The stress tensor will be set to zero',ch10,& 2371 & ' ---' 2372 call wrtout(std_out,message,'COLL') 2373 call wrtout(ab_out,message,'COLL') 2374 effective_potential%strten = zero 2375 else 2376 effective_potential%strten = inp%strten_reference 2377 end if 2378 end if 2379 2380 if(any(abs(effective_potential%strten(:)) >tol16))then 2381 write(message, '(3a)' )ch10,& 2382 & ' Cartesian components of forces (hartree/bohr)',ch10 2383 call wrtout(ab_out,message,'COLL') 2384 call wrtout(std_out, message,'COLL') 2385 do ii = 1, natom 2386 write(message, '(I4,a,3(e16.8))' ) & 2387 & ii,' ',effective_potential%fcart(:,ii) 2388 2389 call wrtout(ab_out,message,'COLL') 2390 call wrtout(std_out, message,'COLL') 2391 end do 2392 2393 write(message, '(a,a)' )ch10,& 2394 & ' Cartesian components of stress tensor (hartree/bohr^3)' 2395 call wrtout(ab_out,message,'COLL') 2396 call wrtout(std_out, message,'COLL') 2397 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2398 & ' sigma(1 1)=',effective_potential%strten(1) / crystal%ucvol,& 2399 & ' sigma(3 2)=',effective_potential%strten(4) / crystal%ucvol 2400 call wrtout(ab_out,message,'COLL') 2401 call wrtout(std_out, message,'COLL') 2402 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2403 & ' sigma(2 2)=',effective_potential%strten(2) / crystal%ucvol,& 2404 & ' sigma(3 1)=',effective_potential%strten(5) / crystal%ucvol 2405 call wrtout(ab_out,message,'COLL') 2406 call wrtout(std_out, message,'COLL') 2407 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2408 & ' sigma(3 3)=',effective_potential%strten(3) / crystal%ucvol,& 2409 & ' sigma(2 1)=',effective_potential%strten(6) / crystal%ucvol 2410 call wrtout(ab_out,message,'COLL') 2411 call wrtout(std_out, message,'COLL') 2412 write(message, '(a)' ) ' ' 2413 call wrtout(ab_out,message,'COLL') 2414 call wrtout(std_out, message,'COLL') 2415 end if 2416 2417 !********************************************************************** 2418 ! Elastic tensors at Gamma Point 2419 !********************************************************************** 2420 write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,& 2421 & ' Extraction of the clamped elastic tensor (unit:10^2GPa)',ch10 2422 call wrtout(std_out,message,'COLL') 2423 call wrtout(ab_out,message,'COLL') 2424 2425 ! look after the blok no.iblok that contains the elastic tensor 2426 qphon(:,1)=zero 2427 qphnrm(1)=zero 2428 rfphon(1:2)=0 2429 rfelfd(1:2)=0 2430 rfstrs(1:2)=3 ! Need uniaxial both stresses and shear stresses 2431 rftyp=1 ! Blocks obtained by a non-stationary formulation. 2432 ! for both diagonal and shear parts 2433 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 2434 2435 if (iblok /=0) then 2436 ! extraction of the elastic constants from the blkvals (GPa) 2437 do ivarA=1,6 2438 do ivarB=1,6 2439 ! because the elastic constant is 6*6, 2440 ! so we should judge if the idir is larger than 3 2441 ! or not 2442 if(ivarA>3) then 2443 idir1=ivarA-3 2444 ipert1=natom+4 !for the shear modulus 2445 else if(ivarA<=3) then 2446 idir1=ivarA 2447 ipert1=natom+3 !for the diagonal part 2448 end if 2449 if(ivarB>3) then 2450 idir2=ivarB-3 2451 ipert2=natom+4 !for the shear modulus 2452 else if(ivarB<=3) then 2453 idir2=ivarB 2454 ipert2=natom+3 !for the diagonal part 2455 end if 2456 elast_clamped(ivarA,ivarB) = blkval(1,idir1,ipert1,idir2,ipert2,iblok) 2457 end do 2458 end do 2459 fact=HaBohr3_GPa / crystal%ucvol 2460 do ivarA=1,6 2461 write(message,'(6f12.7)')elast_clamped(ivarA,1)*fact/100.00_dp,& 2462 & elast_clamped(ivarA,2)*fact/100.00_dp,& 2463 & elast_clamped(ivarA,3)*fact/100.00_dp,& 2464 & elast_clamped(ivarA,4)*fact/100.00_dp,& 2465 & elast_clamped(ivarA,5)*fact/100.00_dp,& 2466 & elast_clamped(ivarA,6)*fact/100.00_dp 2467 call wrtout(std_out,message,'COLL') 2468 call wrtout(ab_out,message,'COLL') 2469 end do 2470 2471 ! Set the clamped tensor into the effective potentiel 2472 effective_potential%harmonics_terms%elastic_constants = elast_clamped 2473 2474 else 2475 2476 write(message,'(3a)')ch10,& 2477 & ' Warning : Elastic Tensor is set to zero (not available in the DDB)' 2478 call wrtout(std_out,message,'COLL') 2479 call wrtout(ab_out,message,'COLL') 2480 2481 ! Set the clamped tensor to zero into the effective potentiel (not available in the DDB) 2482 effective_potential%harmonics_terms%elastic_constants = zero 2483 end if 2484 2485 !********************************************************************** 2486 ! Acoustic Sum Rule 2487 !*************************************************************************** 2488 ! ASR-correction (d2asr) has to be determined here from the Dynamical matrix at Gamma. 2489 ABI_MALLOC(d2asr,(2,3,natom,3,natom)) 2490 2491 write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,& 2492 & ' Calculation of acoustic sum rule',ch10 2493 call wrtout(std_out,message,'COLL') 2494 call wrtout(ab_out,message,'COLL') 2495 2496 ! Find the Gamma block in the DDB (no need for E-field entries) 2497 qphon(:,1)=zero 2498 qphnrm(1)=zero 2499 rfphon(1:2)=1 2500 rfelfd(:)=0 2501 rfstrs(:)=0 2502 rftyp=inp%rfmeth 2503 2504 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 2505 2506 d2asr = zero 2507 if (iblok /=0) then 2508 call asria_calc(inp%asr,d2asr,ddb%val(:,:,iblok),ddb%mpert,ddb%natom) 2509 end if 2510 2511 ! Acoustic Sum Rule 2512 ! In case the interatomic forces are not calculated, the 2513 ! ASR-correction (asrq0%d2asr) has to be determined here from the Dynamical matrix at Gamma. 2514 asrq0 = ddb%get_asrq0(inp%asr, inp%rfmeth, crystal%xcart) 2515 2516 !********************************************************************** 2517 ! Interatomic Forces Calculation 2518 !********************************************************************** 2519 ! ifc to be calculated for interpolation 2520 write(message, '(a,a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,& 2521 & ' Calculation of the interatomic forces from DDB',ch10 2522 call wrtout(std_out,message,'COLL') 2523 call wrtout(ab_out,message,'COLL') 2524 2525 call ifc%init(crystal,ddb,inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,& 2526 & inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,effective_potential%harmonics_terms%zeff,qdrp_cart,& 2527 & inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm) 2528 2529 !*************************************************************************** 2530 ! Interpolation of the dynamical matrix for each qpoint from ifc 2531 !*************************************************************************** 2532 2533 ABI_MALLOC(d2cart,(2,3,mpert,3,mpert)) 2534 ABI_MALLOC(displ,(2*3*natom*3*natom)) 2535 ABI_MALLOC(eigval,(3,natom)) 2536 ABI_MALLOC(eigvec,(2,3,natom,3,natom)) 2537 ABI_MALLOC(phfrq,(3*natom)) 2538 2539 ABI_MALLOC(effective_potential%harmonics_terms%dynmat,(2,3,natom,3,natom,inp%nph1l)) 2540 ABI_MALLOC(effective_potential%harmonics_terms%phfrq,(3*natom,inp%nph1l)) 2541 ABI_MALLOC(effective_potential%harmonics_terms%qpoints,(3,inp%nph1l)) 2542 2543 write(message,'(a,(80a),3a)')ch10,('=',ii=1,80),ch10,ch10,& 2544 & ' Calculation of dynamical matrix for each ph1l points ' 2545 call wrtout(ab_out,message,'COLL') 2546 call wrtout(std_out,message,'COLL') 2547 2548 !Transfer value in effective_potential structure 2549 effective_potential%harmonics_terms%nqpt = inp%nph1l 2550 effective_potential%harmonics_terms%qpoints(:,:) = inp%qph1l(:,:) 2551 2552 ! Store the highest frequency 2553 max_phfq = zero 2554 2555 do iphl1=1,inp%nph1l 2556 2557 ! Initialisation of the phonon wavevector 2558 qphon(:,1)=inp%qph1l(:,iphl1) 2559 if (inp%nph1l /= 0) qphnrm(1) = inp%qnrml1(iphl1) 2560 2561 ! Get d2cart using the interatomic forces and the 2562 ! long-range coulomb interaction through Ewald summation 2563 call gtdyn9(ddb%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip,ifc%dyewq0,d2cart,crystal%gmet,& 2564 & ddb%gprim,mpert,natom,ifc%nrpt,qphnrm(1),qphon(:,1),crystal%rmet,ddb%rprim,ifc%rpt,& 2565 & ifc%trans,crystal%ucvol,ifc%wghatm,crystal%xred,zeff,qdrp_cart,ifc%ewald_option,xmpi_comm_self) 2566 2567 ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix 2568 call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,crystal%indsym,& 2569 & mpert,crystal%nsym,natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),qphon,& 2570 & crystal%rprimd,inp%symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol) 2571 2572 ! Write the phonon frequencies 2573 call dfpt_prtph(displ,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon) 2574 2575 ! Store the highest frequency in cmm-1 2576 max_phfq = max(maxval(phfrq*Ha_cmm1),max_phfq) 2577 2578 effective_potential%harmonics_terms%dynmat(:,:,:,:,:,iphl1) = d2cart(:,:,:natom,:,:natom) 2579 effective_potential%harmonics_terms%phfrq(:,iphl1) = phfrq(:) * Ha_cmm1 2580 2581 end do 2582 2583 write(message, '(2a,f15.7,a)' ) ch10,& 2584 & ' The highest frequency found is ',max_phfq,' cm-1' 2585 call wrtout(std_out,message,'COLL') 2586 2587 ABI_FREE(d2cart) 2588 ABI_FREE(displ) 2589 ABI_FREE(eigval) 2590 ABI_FREE(eigvec) 2591 ABI_FREE(phfrq) 2592 2593 !********************************************************************** 2594 ! Transfert inter-atomic forces constants in reduced coordinates 2595 !********************************************************************** 2596 2597 !Reorder cell from canonical coordinates to reduced coordinates (for multibinit) 2598 !store the number of ifc before rearrangement 2599 2600 ! Store the sum of the weight of IFC for the final check 2601 wcount1 = 0 2602 do irpt=1,ifc%nrpt 2603 wcount1 = wcount1 + sum(ifc%wghatm(:,:,irpt)) 2604 end do 2605 2606 !Set the maximum and the miminum for the bound of the cell 2607 max1 = maxval(ifc%cell(1,:)); min1 = minval(ifc%cell(1,:)) 2608 max2 = maxval(ifc%cell(2,:)); min2 = minval(ifc%cell(2,:)) 2609 max3 = maxval(ifc%cell(3,:)); min3 = minval(ifc%cell(3,:)) 2610 cell_number(1) = max1 - min1 + 1 2611 cell_number(2) = max2 - min2 + 1 2612 cell_number(3) = max3 - min3 + 1 2613 2614 ! set the new number of cell, sometimes, in canonical coordinates, 2615 ! some cell are delete but they exist in reduced coordinates. 2616 nrpt_new = product(cell_number(:)) 2617 2618 ! Allocate temporary array 2619 ABI_MALLOC(atmfrc_red,(3,natom,3,natom,nrpt_new)) 2620 ABI_MALLOC(wghatm_red,(natom,natom,nrpt_new)) 2621 ABI_MALLOC(cell_red,(3,nrpt_new)) 2622 2623 wghatm_red(:,:,:) = zero 2624 2625 if(iam_master)then 2626 do ia=1,natom 2627 do ib=1,natom 2628 2629 ! Simple Lattice 2630 if (inp%brav==1) then 2631 ! In this case, it is better to work in reduced coordinates 2632 ! As rcan is in canonical coordinates, => multiplication by gprim 2633 do ii=1,3 2634 red(1,ii)= ifc%rcan(1,ia)*ddb%gprim(1,ii) + & 2635 & ifc%rcan(2,ia)*ddb%gprim(2,ii) + & 2636 & ifc%rcan(3,ia)*ddb%gprim(3,ii) 2637 red(2,ii)= ifc%rcan(1,ib)*ddb%gprim(1,ii) + & 2638 ifc%rcan(2,ib)*ddb%gprim(2,ii) + & 2639 & ifc%rcan(3,ib)*ddb%gprim(3,ii) 2640 end do 2641 end if 2642 2643 ! Get the shift of cell 2644 shift(:) = int(anint(red(2,:) - crystal%xred(:,ib)) - anint(red(1,:) - crystal%xred(:,ia))) 2645 2646 do irpt=1,ifc%nrpt 2647 2648 cell2(:)= int(ifc%cell(:,irpt) + shift(:)) 2649 2650 ! Use boundary condition to get the right cell 2651 if (cell2(1) < min1 .and. cell2(1) < max1) then 2652 cell2(1) = cell2(1) + cell_number(1) 2653 else if (cell2(1) > min1 .and. cell2(1) > max1) then 2654 cell2(1) = cell2(1) - cell_number(1) 2655 end if 2656 2657 if (cell2(2) < min2 .and. cell2(2) < max2) then 2658 cell2(2) = cell2(2) + cell_number(2) 2659 else if (cell2(2) > min2 .and. cell2(2) > max2) then 2660 cell2(2) = cell2(2) - cell_number(2) 2661 end if 2662 2663 if (cell2(3) < min3 .and. cell2(3) < max3) then 2664 cell2(3) = cell2(3) + cell_number(3) 2665 else if (cell2(3) > min3 .and. cell2(3) > max3) then 2666 cell2(3) = cell2(3) - cell_number(3) 2667 end if 2668 2669 irpt2=1 2670 do i1=min1,max1 2671 do i2=min2,max2 2672 do i3=min3,max3 2673 if (i1 == cell2(1) .and.& 2674 i2 == cell2(2) .and.& 2675 i3 == cell2(3)) then 2676 wghatm_red(ia,ib,irpt2) = ifc%wghatm(ia,ib,irpt) 2677 atmfrc_red(:,ia,:,ib,irpt2) = ifc%atmfrc(:,ia,:,ib,irpt) 2678 cell_red(1,irpt2) = i1 2679 cell_red(2,irpt2) = i2 2680 cell_red(3,irpt2) = i3 2681 end if 2682 irpt2 = irpt2 + 1 2683 end do 2684 end do 2685 end do 2686 end do 2687 end do 2688 end do 2689 end if 2690 2691 call xmpi_bcast(atmfrc_red,master, comm, ierr) 2692 call xmpi_bcast(wghatm_red,master, comm, ierr) 2693 call xmpi_bcast(cell_red,master, comm, ierr) 2694 2695 ! Copy ifc into effective potential 2696 ! !!Warning eff_pot%ifcs only contains atmfrc,short_atmfrc,ewald_atmfrc,,nrpt and cell!! 2697 ! rcan,ifc%rpt,wghatm and other quantities 2698 ! are not needed for effective potential!!! 2699 call ifc%free() 2700 call effective_potential%harmonics_terms%ifcs%free() 2701 2702 ! Only conserve the necessary points in rpt 2703 nrpt_new2 = 0 2704 do irpt = 1, nrpt_new 2705 if (abs(sum(wghatm_red(:,:,irpt))) >tol16) then 2706 nrpt_new2 = nrpt_new2 + 1 2707 end if 2708 end do 2709 2710 ! Set the new number of rpt 2711 effective_potential%harmonics_terms%ifcs%nrpt = nrpt_new2 2712 2713 ! Allocation of the final arrays 2714 ABI_MALLOC(effective_potential%harmonics_terms%ifcs%atmfrc,(3,natom,3,natom,nrpt_new2)) 2715 ABI_MALLOC(effective_potential%harmonics_terms%ifcs%short_atmfrc,(3,natom,3,natom,nrpt_new2)) 2716 ABI_MALLOC(effective_potential%harmonics_terms%ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt_new2)) 2717 ABI_MALLOC(effective_potential%harmonics_terms%ifcs%cell,(3,nrpt_new2)) 2718 ABI_MALLOC(effective_potential%harmonics_terms%ifcs%wghatm,(natom,natom,nrpt_new2)) 2719 2720 irpt2 = 0 2721 do irpt = 1,nrpt_new 2722 if (abs(sum(wghatm_red(:,:,irpt))) > tol16) then 2723 irpt2 = irpt2 + 1 2724 ! Apply weight on each R point 2725 do ia=1,effective_potential%crystal%natom 2726 do ib=1,effective_potential%crystal%natom 2727 atmfrc_red(:,ia,:,ib,irpt) = atmfrc_red(:,ia,:,ib,irpt)*wghatm_red(ia,ib,irpt) 2728 end do 2729 end do 2730 effective_potential%harmonics_terms%ifcs%cell(:,irpt2) = cell_red(:,irpt) 2731 effective_potential%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt2) = atmfrc_red(:,:,:,:,irpt) 2732 if (inp%dipdip == 1) then 2733 effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)=& 2734 & atmfrc_red(:,:,:,:,irpt) 2735 else 2736 effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2) = zero 2737 end if 2738 effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)=atmfrc_red(:,:,:,:,irpt) 2739 effective_potential%harmonics_terms%ifcs%ewald_atmfrc(:,:,:,:,irpt2) = zero 2740 effective_potential%harmonics_terms%ifcs%wghatm(:,:,irpt2) = wghatm_red(:,:,irpt) 2741 end if 2742 end do 2743 2744 2745 ABI_FREE(atmfrc_red) 2746 ABI_FREE(wghatm_red) 2747 ABI_FREE(cell_red) 2748 2749 ! Final check 2750 wcount2 = 0 2751 do irpt = 1, effective_potential%harmonics_terms%ifcs%nrpt 2752 wcount2 = wcount2 + sum(effective_potential%harmonics_terms%ifcs%wghatm(:,:,irpt)) 2753 end do 2754 2755 if (abs(wcount1-wcount2)/(wcount1+wcount2)>tol8) then 2756 write(message,'(2a,es15.4,a,es15.4,a,es15.4)')'The total wghatm has changed',ch10,& 2757 & wcount1,' before and ', wcount2, ' now, difference being ',wcount1-wcount2 2758 ABI_BUG(message) 2759 end if 2760 2761 2762 !********************************************************************** 2763 ! Internal strain tensors at Gamma point 2764 !********************************************************************** 2765 write(message, '(a,a,(80a),a,a,a)') ch10,('=',ii=1,80),ch10,ch10,& 2766 & ' Calculation of the internal-strain tensor' 2767 call wrtout(std_out,message,'COLL') 2768 call wrtout(ab_out,message,'COLL') 2769 ABI_MALLOC(instrain,(3*natom,6)) 2770 ! looking after the no. of blok that contains the internal strain tensor 2771 qphon(:,1)=zero 2772 qphnrm(1)=zero 2773 rfphon(1:2)=0 2774 rfelfd(1:2)=0 2775 rfstrs(1:2)=3 2776 rftyp=1 2777 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 2778 2779 ABI_MALLOC(effective_potential%harmonics_terms%strain_coupling,(6,3,natom)) 2780 effective_potential%harmonics_terms%strain_coupling = zero 2781 2782 if (iblok /=0) then 2783 2784 ! then print the internal strain tensor (only the force one) 2785 prt_internalstr=1 2786 call ddb_internalstr(inp%asr,ddb%val,d2asr,iblok,instrain,& 2787 & ab_out,mpert,natom,nblok,prt_internalstr) 2788 2789 do ipert1=1,6 2790 do ipert2=1,natom 2791 do idir2=1,3 2792 ii=3*(ipert2-1)+idir2 2793 effective_potential%harmonics_terms%strain_coupling(ipert1,idir2,ipert2)=& 2794 & (-1.0_dp)*instrain(ii,ipert1) 2795 end do 2796 end do 2797 end do 2798 else 2799 write(message,'(3a)')ch10,& 2800 & ' Warning : Internal strain is set to zero (not available in the DDB)' 2801 call wrtout(std_out,message,'COLL') 2802 call wrtout(ab_out,message,'COLL') 2803 end if 2804 !------------------------------------------------------------------------------------- 2805 ! DEALLOCATION OF ARRAYS 2806 ABI_FREE(blkval) 2807 ABI_FREE(zeff) 2808 ABI_FREE(qdrp_cart) 2809 ABI_FREE(instrain) 2810 ABI_FREE(d2asr) 2811 call asrq0%free() 2812 2813 write(message,'(a)')ch10 2814 call wrtout(std_out,message,'COLL') 2815 call wrtout(ab_out,message,'COLL') 2816 2817 end subroutine system_ddb2effpot
m_effective_potential_file/system_getDimFromXML [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
system_getDimFromXML
FUNCTION
Open xml file of effective potentiel, then reads the variables that must be known in order to dimension the arrays before complete reading
INPUTS
character(len=*) filnam: name of input or output file
OUTPUT
natom=number of atoms ntypat=number of atom types nrpt =number of real space points used to integrate IFC
SOURCE
1050 subroutine system_getDimFromXML(filename,natom,ntypat,nph1l,nrpt) 1051 1052 !Arguments ------------------------------------ 1053 !scalars 1054 character(len=fnlen),intent(in) :: filename 1055 integer, intent(out) :: natom,ntypat,nph1l,nrpt 1056 !arrays 1057 !Local variables------------------------------- 1058 !scalar 1059 integer :: nrpt1,nrpt2 1060 real :: itypat 1061 character(len=500) :: message 1062 #ifndef HAVE_XML 1063 integer :: funit = 1,ios = 0 1064 integer :: iatom 1065 logical :: found 1066 character (len=XML_RECL) :: line,readline 1067 character (len=XML_RECL) :: strg,strg1 1068 #endif 1069 !arrays 1070 #ifndef HAVE_XML 1071 integer,allocatable :: typat(:) 1072 #endif 1073 1074 ! ************************************************************************* 1075 1076 !Open the atomicdata XML file for reading 1077 write(message,'(5a)') ' system_getDimFromXML :',& 1078 & '-Opening the file ',trim(filename),' to read dimensions',& 1079 & ' (before initialisation)' 1080 1081 call wrtout(std_out,message,'COLL') 1082 1083 natom = 0 1084 ntypat= 0 1085 nph1l = 0 1086 nrpt = 0 1087 nrpt1 = 0 1088 nrpt2 = 0 1089 itypat= 0 1090 1091 !Open the atomicdata XML file for reading 1092 1093 #if defined HAVE_XML 1094 !Read with libxml 1095 call effpot_xml_getDimSystem(char_f2c(trim(filename)),natom,ntypat,nph1l,nrpt1,nrpt2) 1096 #else 1097 !Read by hand 1098 1099 !Start a reading loop 1100 found=.false. 1101 1102 if (open_file(filename,message,unit=funit,form="formatted",status="old",& 1103 & action="read") /= 0) then 1104 ABI_ERROR(message) 1105 end if 1106 1107 !First parse to know the number of atoms 1108 do while ((ios==0).and.(.not.found)) 1109 read(funit,'(a)',iostat=ios) readline 1110 if(ios ==0)then 1111 call rmtabfromline(readline) 1112 line=adjustl(readline) 1113 1114 !Need test with char(9) because the old version of XML file 1115 !from old script includes tarbulation at the begining of each line 1116 if (line(1:5)==char(60)//'atom') then 1117 natom=natom+1 1118 cycle 1119 end if 1120 1121 if (line(1:21)==char(60)//'local_force_constant') then 1122 nrpt1 = nrpt1+1 1123 cycle 1124 end if 1125 1126 if (line(1:21)==char(60)//'total_force_constant') then 1127 nrpt2 = nrpt2+1 1128 cycle 1129 end if 1130 1131 if (line(1:7)==char(60)//'qpoint') then 1132 nph1l = nph1l+1 1133 cycle 1134 end if 1135 end if 1136 end do 1137 1138 !second parse to get the number of typat 1139 ABI_MALLOC(typat,(natom)) 1140 typat = 0 1141 iatom = 0 1142 1143 rewind(funit) 1144 !Start a reading loop 1145 ios = 0 1146 found = .false. 1147 1148 do while ((ios==0).and.(.not.found)) 1149 read(funit,'(a)',iostat=ios) readline 1150 if(ios == 0)then 1151 call rmtabfromline(readline) 1152 line=adjustl(readline) 1153 1154 if (line(1:5)==char(60)//'atom') then 1155 iatom = iatom + 1 1156 call rdfromline("mass",line,strg) 1157 strg1=trim(strg) 1158 read(unit=strg1,fmt=*) itypat 1159 if (.not.any(typat==int(itypat))) then 1160 ntypat= ntypat+1 1161 end if 1162 typat(iatom) = int(itypat) 1163 end if 1164 1165 if (line(1:6)==char(60)//'local') then 1166 found=.true. 1167 end if 1168 end if 1169 end do 1170 1171 close(funit) 1172 ABI_FREE(typat) 1173 1174 #endif 1175 1176 !Check the RPT 1177 if (nrpt2/=nrpt1) then 1178 if(nrpt1> 0 .and. nrpt2== 0) then 1179 continue; 1180 else if (nrpt1==0.and.nrpt2>=0) then 1181 write(message, '(5a)' )ch10,& 1182 & ' WARNING: the number of local IFC is set to 0 ',ch10,& 1183 & ' Dipdip must be set to zero',ch10 1184 call wrtout(std_out,message,'COLL') 1185 else if (nrpt2 > nrpt1) then 1186 write(message, '(2a,I0,3a,I0,5a)' )ch10,& 1187 & ' WARNING: the number of total IFC (',nrpt2,') is not equal to the ',ch10,& 1188 & ' the number of short range IFC (',nrpt1,') in ',filename,ch10,& 1189 & ' the missing ifc will be set to zero',ch10 1190 call wrtout(std_out,message,'COLL') 1191 else if(nrpt1>nrpt2)then 1192 write(message, '(2a,I0,3a,I0,5a)' )ch10,& 1193 & ' The number of total IFC (',nrpt2,') is inferior to ',ch10,& 1194 & ' the number of short range IFC (',nrpt1,') in ',filename,ch10,& 1195 & ' This is not possible',ch10 1196 ABI_BUG(message) 1197 end if 1198 end if 1199 1200 !nrpt is the max between local and total: 1201 nrpt = max(nrpt1,nrpt2) 1202 1203 1204 end subroutine system_getDimFromXML
m_effective_potential_file/system_xml2effpot [ Functions ]
[ Top ] [ m_effective_potential_file ] [ Functions ]
NAME
system_xml2effpot
FUNCTION
Open xml file of effective potentiel, then reads the variables and store them in effective potentential type
INPUTS
eff_pot<type(effective_potential_type)> = datatype with all the information for effective potential comm=MPI communicator character(len=*) filnam: name of input or output file strcpling = optional,logical to disable the strcpling
OUTPUT
eff_pot<type(effective_potential_type)> = datatype with all the information for effective potential
SOURCE
1226 subroutine system_xml2effpot(eff_pot,filename,comm,strcpling) 1227 1228 use m_atomdata 1229 use m_multibinit_dataset, only : multibinit_dtset_type 1230 use m_ab7_symmetry 1231 1232 !Arguments ------------------------------------ 1233 !scalars 1234 character(len=*),intent(in) :: filename 1235 integer, intent(in) :: comm 1236 integer, optional,intent(in) :: strcpling 1237 !arrays 1238 type(effective_potential_type), intent(inout) :: eff_pot 1239 1240 !Local variables------------------------------- 1241 !scalar 1242 integer :: ierr,ii,itypat,my_rank,msym,natom,ncoeff,nrpt,nrpt_scoupling 1243 integer :: ntypat,nph1l,nptsym,npsp,nproc,nsym,space_group,timrev,use_inversion,voigt 1244 real(dp):: energy,tolsym,ucvol 1245 character(len=500) :: message 1246 integer,parameter :: master=0 1247 logical :: has_anharmonics 1248 logical :: iam_master 1249 #ifndef HAVE_XML 1250 integer :: funit = 1,ios=0 1251 integer :: iatom,iamu,iph1l,irpt,irpt1,irpt2,irpt3,jj,mu,nu 1252 real(dp):: amu 1253 logical :: found,found2,short_range,total_range 1254 character (len=XML_RECL) :: line,readline 1255 character (len=XML_RECL) :: strg,strg1 1256 logical :: has_straincoupling 1257 #endif 1258 !arrays 1259 integer :: bravais(11) 1260 integer,allocatable :: typat(:) 1261 integer,allocatable :: symrel(:,:,:),symafm(:),ptsymrel(:,:,:) 1262 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 1263 real(dp) :: elastic_constants(6,6),elastic3rd(6,6,6),epsilon_inf(3,3) 1264 real(dp),allocatable :: all_amu(:),cell_local(:,:),cell_total(:,:) 1265 real(dp),allocatable :: elastic_displacement(:,:,:,:),dynmat(:,:,:,:,:,:) 1266 real(dp),allocatable :: local_atmfrc(:,:,:,:,:),total_atmfrc(:,:,:,:,:) 1267 real(dp),allocatable :: spinat(:,:),strain_coupling(:,:,:),phfrq(:,:),qph1l(:,:),tnons(:,:) 1268 real(dp),allocatable :: xcart(:,:),xred(:,:),zeff(:,:,:),znucl(:),zion(:) 1269 character(len=132),allocatable :: title(:) 1270 type(ifc_type) :: ifcs 1271 type(ifc_type),dimension(:),allocatable :: phonon_strain 1272 type(crystal_t) :: crystal 1273 type(atomdata_t) :: atom 1274 #ifdef HAVE_XML 1275 real(dp),allocatable :: phonon_strain_atmfrc(:,:,:,:,:) 1276 integer,allocatable :: phonon_straincell(:,:) 1277 #endif 1278 #ifndef HAVE_XML 1279 real(dp),allocatable :: work2(:,:) 1280 #endif 1281 1282 ! ************************************************************************* 1283 1284 1285 !Open the atomicdata XML file for reading 1286 write(message,'(a,a)')'-Opening the file ',filename 1287 1288 call wrtout(ab_out,message,'COLL') 1289 call wrtout(std_out,message,'COLL') 1290 1291 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1292 iam_master = (my_rank == master) 1293 1294 !Get Dimention of system and allocation/initialisation of array 1295 call effective_potential_file_getDimSystem(filename,comm,natom,ntypat,nph1l,nrpt) 1296 gmet= zero; gprimd = zero; rmet = zero; rprimd = zero 1297 elastic_constants = zero; epsilon_inf = zero; ncoeff = 0 1298 ABI_MALLOC(all_amu,(ntypat)) 1299 ABI_MALLOC(cell_local,(3,nrpt)) 1300 ABI_MALLOC(cell_total,(3,nrpt)) 1301 ABI_MALLOC(elastic_displacement,(6,6,3,natom)) 1302 ABI_MALLOC(ifcs%atmfrc,(3,natom,3,natom,nrpt)) 1303 ABI_MALLOC(ifcs%cell,(3,nrpt)) 1304 ABI_MALLOC(ifcs%short_atmfrc,(3,natom,3,natom,nrpt)) 1305 ABI_MALLOC(ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt)) 1306 ABI_MALLOC(strain_coupling,(6,3,natom)) 1307 ABI_MALLOC(total_atmfrc,(3,natom,3,natom,nrpt)) 1308 ABI_MALLOC(local_atmfrc,(3,natom,3,natom,nrpt)) 1309 ABI_MALLOC(dynmat,(2,3,natom,3,natom,nph1l)) 1310 ABI_MALLOC(typat,(natom)) 1311 ABI_MALLOC(phfrq,(3*natom,nph1l)) 1312 ABI_MALLOC(qph1l,(3,nph1l)) 1313 ABI_MALLOC(xcart,(3,natom)) 1314 ABI_MALLOC(xred,(3,natom)) 1315 ABI_MALLOC(zeff,(3,3,natom)) 1316 ABI_MALLOC(zion,(ntypat)) 1317 ABI_MALLOC(znucl,(ntypat)) 1318 1319 ABI_MALLOC(phonon_strain,(6)) 1320 nrpt_scoupling = 0 1321 do ii = 1,6 1322 ! Get The size of the strainPhonon-coupling 1323 call effective_potential_file_getDimStrainCoupling(filename,nrpt_scoupling,ii-1) 1324 ABI_MALLOC(phonon_strain(ii)%atmfrc,(3,natom,3,natom,nrpt_scoupling)) 1325 ABI_MALLOC(phonon_strain(ii)%cell,(3,nrpt_scoupling)) 1326 phonon_strain(ii)%nrpt = nrpt_scoupling 1327 phonon_strain(ii)%atmfrc = zero 1328 phonon_strain(ii)%cell = 0 1329 end do 1330 1331 all_amu(:) = zero 1332 dynmat(:,:,:,:,:,:) = zero 1333 cell_local(:,:) = 99D99 1334 cell_total(:,:) = 99D99 1335 elastic3rd(:,:,:) = zero 1336 elastic_displacement(:,:,:,:) = zero 1337 ifcs%nrpt = nrpt 1338 ifcs%atmfrc(:,:,:,:,:) = zero 1339 ifcs%cell(:,:) = 0 1340 ifcs%ewald_atmfrc(:,:,:,:,:) = zero 1341 ifcs%short_atmfrc(:,:,:,:,:) = zero 1342 strain_coupling(:,:,:) = zero 1343 phfrq = zero 1344 qph1l = 0 1345 xcart = zero 1346 zeff = zero 1347 znucl = zero 1348 1349 if(iam_master)then 1350 !Open the atomicdata XML file for reading 1351 #if defined HAVE_XML 1352 1353 write(message,'(a,a,a,a)')'-Reading the file ',trim(filename),& 1354 & ' with LibXML library' 1355 1356 call wrtout(ab_out,message,'COLL') 1357 call wrtout(std_out,message,'COLL') 1358 1359 ! Read with libxml library 1360 call effpot_xml_readSystem(char_f2c(trim(filename)),natom,ntypat,nrpt,nph1l,all_amu,& 1361 & ifcs%atmfrc,ifcs%cell,dynmat,elastic_constants,energy,& 1362 & epsilon_inf,ifcs%ewald_atmfrc,phfrq,rprimd,qph1l,& 1363 & ifcs%short_atmfrc,typat,xcart,zeff) 1364 1365 ! convert atomic mass unit to znucl 1366 do itypat=1,ntypat 1367 do ii=1,103 1368 call atomdata_from_znucl(atom,real(ii,dp)) 1369 if (abs((real(atom%amu,sp)-real(all_amu(itypat),sp))& 1370 & /real(all_amu(itypat),sp)*100)<0.1) then 1371 znucl(itypat) = atom%znucl 1372 exit 1373 end if 1374 end do 1375 end do 1376 1377 ! Get the Phonon Strain coupling 1378 do voigt = 1,6 1379 nrpt_scoupling = phonon_strain(voigt)%nrpt 1380 ABI_MALLOC(phonon_straincell,(3,nrpt_scoupling)) 1381 ABI_MALLOC(phonon_strain_atmfrc,(3,natom,3,natom,nrpt_scoupling)) 1382 1383 ! Get The value 1384 call effpot_xml_readStrainCoupling(char_f2c(trim(filename)),natom,nrpt_scoupling,(voigt-1),& 1385 & elastic3rd(voigt,:,:),elastic_displacement(voigt,:,:,:),& 1386 & strain_coupling(voigt,:,:),& 1387 & phonon_strain_atmfrc,phonon_straincell) 1388 1389 ! Check if the 3rd order strain_coupling is present 1390 has_anharmonics = .FALSE. 1391 if(any(elastic3rd>tol10).or.any(elastic_displacement>tol10)) has_anharmonics = .TRUE. 1392 phonon_strain(voigt)%atmfrc(:,:,:,:,:) = phonon_strain_atmfrc(:,:,:,:,:) 1393 phonon_strain(voigt)%cell(:,:) = phonon_straincell(:,:) 1394 if(any(phonon_strain(voigt)%atmfrc > tol10)) has_anharmonics = .TRUE. 1395 1396 ABI_FREE(phonon_straincell) 1397 ABI_FREE(phonon_strain_atmfrc) 1398 end do 1399 #else 1400 1401 ! Read by hand 1402 write(message,'(a,a,a,a)')'-Reading the file ',trim(filename),& 1403 & ' with Fortran' 1404 1405 call wrtout(ab_out,message,'COLL') 1406 call wrtout(std_out,message,'COLL') 1407 1408 if (open_file(filename,message,unit=funit,form="formatted",& 1409 & status="old",action="read") /= 0) then 1410 ABI_ERROR(message) 1411 end if 1412 1413 !Start a reading loop in fortran 1414 rewind(unit=funit) 1415 found=.false. 1416 1417 iatom = 1 1418 iamu = 1 1419 itypat = 1 1420 irpt = 1 1421 irpt1 = 0 1422 irpt2 = 0 1423 iph1l = 1 1424 amu = zero 1425 short_range = .false. 1426 total_range = .false. 1427 has_straincoupling = .FALSE. 1428 do while ((ios==0).and.(.not.found)) 1429 read(funit,'(a)',iostat=ios) readline 1430 if(ios == 0)then 1431 call rmtabfromline(readline) 1432 line=adjustl(readline) 1433 if (.not.has_straincoupling) then 1434 1435 if ((line(1:7)=='<energy')) then 1436 call rdfromline_value('energy',line,strg) 1437 if (strg/="") then 1438 strg1=trim(strg) 1439 read(strg1,*) energy 1440 else 1441 read(funit,'(a)',iostat=ios) readline 1442 call rmtabfromline(readline) 1443 line=adjustl(readline) 1444 call rdfromline_value('energy',line,strg) 1445 if (strg/="") then 1446 strg1=trim(strg) 1447 else 1448 strg1=trim(line) 1449 end if 1450 read(strg1,*) energy 1451 end if 1452 cycle 1453 end if 1454 1455 if ((line(1:10)=='<unit_cell')) then 1456 call rdfromline_value('unit_cell',line,strg) 1457 if (strg/="") then 1458 strg1=trim(strg) 1459 read(strg1,*) (rprimd(1,mu),mu=1,3) 1460 read(funit,*) (rprimd(2,mu),mu=1,3) 1461 else 1462 do nu=1,2 1463 read(funit,*) (rprimd(nu,mu),mu=1,3) 1464 end do 1465 end if 1466 read(funit,'(a)',iostat=ios) readline 1467 call rmtabfromline(readline) 1468 line=adjustl(readline) 1469 call rdfromline_value('unit_cell',line,strg) 1470 if (strg/="") then 1471 strg1=trim(strg) 1472 else 1473 strg1=trim(line) 1474 end if 1475 read(strg1,*) (rprimd(3,mu),mu=1,3) 1476 cycle 1477 end if 1478 1479 if ((line(1:12)=='<epsilon_inf')) then 1480 call rdfromline_value('epsilon_inf',line,strg) 1481 if (strg/="") then 1482 strg1=trim(strg) 1483 read(strg1,*) (epsilon_inf(mu,1),mu=1,3) 1484 read(funit,*) (epsilon_inf(mu,2),mu=1,3) 1485 else 1486 do nu=1,2 1487 read(funit,*) (epsilon_inf(mu,nu),mu=1,3) 1488 end do 1489 end if 1490 read(funit,'(a)',iostat=ios) readline 1491 call rmtabfromline(readline) 1492 line=adjustl(readline) 1493 call rdfromline_value('epsilon_inf',line,strg) 1494 if (strg/="") then 1495 strg1=trim(strg) 1496 else 1497 strg1=trim(line) 1498 end if 1499 read(strg1,*) (epsilon_inf(mu,3),mu=1,3) 1500 cycle 1501 end if 1502 1503 if ((line(1:8)=='<elastic')) then 1504 call rdfromline_value('elastic',line,strg) 1505 if (strg/="") then 1506 strg1=trim(strg) 1507 read(strg1,*) (elastic_constants(mu,1),mu=1,6) 1508 do nu=2,5 1509 read(funit,*) (elastic_constants(mu,nu),mu=1,6) 1510 end do 1511 else 1512 do nu=1,5 1513 read(funit,*) (elastic_constants(mu,nu),mu=1,6) 1514 end do 1515 end if 1516 read(funit,'(a)',iostat=ios) readline 1517 call rmtabfromline(readline) 1518 line=adjustl(readline) 1519 call rdfromline_value('elastic',line,strg) 1520 if (strg/="") then 1521 strg1=trim(strg) 1522 else 1523 strg1=trim(line) 1524 end if 1525 read(strg1,*) (elastic_constants(mu,6),mu=1,6) 1526 cycle 1527 end if 1528 1529 if ((line(1:5)=='<atom')) then 1530 call rdfromline("mass",line,strg) 1531 strg1=trim(strg) 1532 read(strg1,*) amu 1533 if (.not.any(abs(all_amu-amu)<tol16)) then 1534 all_amu(iamu) = amu 1535 typat(iatom) = int(amu) 1536 !convert atomic mass unit to znucl 1537 do ii=1,103 1538 call atomdata_from_znucl(atom,real(ii,dp)) 1539 if (abs((real(atom%amu,sp)-real(amu,sp))& 1540 & /real(amu,sp)*100)<0.1) then 1541 znucl(iamu) = atom%znucl 1542 exit 1543 end if 1544 end do 1545 iamu = iamu +1 1546 end if 1547 do itypat=1,ntypat 1548 if(abs(amu-all_amu(itypat))<tol16) then 1549 typat(iatom) = itypat 1550 end if 1551 end do 1552 cycle 1553 end if 1554 1555 if ((line(1:9)=='<position')) then 1556 call rdfromline_value('position',line,strg) 1557 if (strg/="") then 1558 strg1=trim(strg) 1559 read(strg1,*)(xcart(mu,iatom),mu=1,3) 1560 else 1561 read(funit,'(a)',iostat=ios) readline 1562 call rmtabfromline(readline) 1563 line=adjustl(readline) 1564 call rdfromline_value('position',line,strg) 1565 if (strg/="") then 1566 strg1=trim(strg) 1567 else 1568 strg1=trim(line) 1569 end if 1570 read(strg1,*)(xcart(mu,iatom),mu=1,3) 1571 end if 1572 cycle 1573 end if 1574 1575 if ((line(1:11)=='<borncharge')) then 1576 call rdfromline_value('borncharge',line,strg) 1577 if (strg/="") then 1578 strg1=trim(strg) 1579 read(strg1,*) (zeff(mu,1,iatom),mu=1,3) 1580 read(funit,*) (zeff(mu,2,iatom),mu=1,3) 1581 else 1582 do nu=1,2 1583 read(funit,*) (zeff(mu,nu,iatom),mu=1,3) 1584 end do 1585 end if 1586 read(funit,'(a)',iostat=ios) readline 1587 line=adjustl(readline) 1588 call rdfromline_value('borncharge',line,strg) 1589 if (strg/="") then 1590 strg1=trim(strg) 1591 else 1592 strg1=trim(line) 1593 end if 1594 read(strg1,*) (zeff(mu,3,iatom),mu=1,3) 1595 cycle 1596 end if 1597 1598 if ((line(1:7)==char(60)//char(47)//'atom'//char(62))) then 1599 iatom=iatom+1 1600 cycle 1601 end if 1602 1603 if ((line(1:12)=='<local_force')) then 1604 found2 = .False. 1605 irpt1 = irpt1 + 1 1606 do while (.not.found2) 1607 read(funit,'(a)',iostat=ios) readline 1608 call rmtabfromline(readline) 1609 line=adjustl(readline) 1610 if ((line(1:5)=='<data')) then 1611 call rdfromline_value('data',line,strg) 1612 if (strg/="") then 1613 ABI_MALLOC(work2,(3*natom,3*natom)) 1614 strg1=trim(strg) 1615 read(strg1,*) (work2(1,nu),nu=1,3*natom) 1616 do mu=2,3*natom-1 1617 read(funit,*)(work2(mu,nu),nu=1,3*natom) 1618 end do 1619 read(funit,'(a)',iostat=ios) readline 1620 call rmtabfromline(readline) 1621 line=adjustl(readline) 1622 call rdfromline_value('data',line,strg) 1623 if (strg/="") then 1624 strg1=trim(strg) 1625 else 1626 strg1=trim(line) 1627 end if 1628 read(strg1,*) (work2(3*natom,nu),nu=1,3*natom) 1629 local_atmfrc(:,:,:,:,irpt1) = reshape(work2,(/3,natom,3,natom/)) 1630 ABI_FREE(work2) 1631 else 1632 ABI_MALLOC(work2,(3*natom,3*natom)) 1633 do mu=1,3*natom 1634 read(funit,*)(work2(mu,nu),nu=1,3*natom) 1635 end do 1636 local_atmfrc(:,:,:,:,irpt1) = reshape(work2,(/3,natom,3,natom/)) 1637 ABI_FREE(work2) 1638 end if 1639 end if 1640 if ((line(1:5)=='<cell')) then 1641 call rdfromline_value('cell',line,strg) 1642 if (strg/="") then 1643 strg1=trim(strg) 1644 read(strg1,*)(cell_local(mu,irpt1),mu=1,3) 1645 else 1646 read(funit,*)(cell_local(mu,irpt1),mu=1,3) 1647 end if 1648 found2 = .TRUE. 1649 cycle 1650 end if 1651 end do 1652 end if 1653 1654 if ((line(1:12)=='<total_force')) then 1655 irpt2 = irpt2 + 1 1656 found2 = .False. 1657 do while (.not.found2) 1658 read(funit,'(a)',iostat=ios) readline 1659 call rmtabfromline(readline) 1660 line=adjustl(readline) 1661 if ((line(1:5)=='<data')) then 1662 call rdfromline_value('data',line,strg) 1663 if (strg/="") then 1664 ABI_MALLOC(work2,(3*natom,3*natom)) 1665 strg1=trim(strg) 1666 read(strg1,*) (work2(1,nu),nu=1,3*natom) 1667 do mu=2,3*natom-1 1668 read(funit,*)(work2(mu,nu),nu=1,3*natom) 1669 end do 1670 read(funit,'(a)',iostat=ios) readline 1671 call rmtabfromline(readline) 1672 line=adjustl(readline) 1673 call rdfromline_value('data',line,strg) 1674 if (strg/="") then 1675 strg1=trim(strg) 1676 else 1677 strg1=trim(line) 1678 end if 1679 read(strg1,*) (work2(3*natom,nu),nu=1,3*natom) 1680 total_atmfrc(:,:,:,:,irpt2) = reshape(work2,(/3,natom,3,natom/)) 1681 ABI_FREE(work2) 1682 else 1683 ABI_MALLOC(work2,(3*natom,3*natom)) 1684 do mu=1,3*natom 1685 read(funit,*)(work2(mu,nu),nu=1,3*natom) 1686 end do 1687 total_atmfrc(:,:,:,:,irpt2) = reshape(work2,(/3,natom,3,natom/)) 1688 ABI_FREE(work2) 1689 end if 1690 end if 1691 if ((line(1:5)=='<cell')) then 1692 call rdfromline_value('cell',line,strg) 1693 if (strg/="") then 1694 strg1=trim(strg) 1695 read(strg1,*)(cell_total(mu,irpt2),mu=1,3) 1696 else 1697 read(funit,*)(cell_total(mu,irpt2),mu=1,3) 1698 end if 1699 found2 = .TRUE. 1700 cycle 1701 end if 1702 end do 1703 end if 1704 1705 if ((line(1:7)=='<qpoint')) then 1706 call rdfromline_value('qpoint',line,strg) 1707 if (strg/="") then 1708 strg1=trim(strg) 1709 read(strg1,*)(qph1l(mu,iph1l),mu=1,3) 1710 else 1711 read(funit,*) (qph1l(mu,iph1l),mu=1,3) 1712 end if 1713 end if 1714 1715 if ((line(1:12)=='<frequencies')) then 1716 call rdfromline_value('frequencies',line,strg) 1717 if (strg/="") then 1718 strg1=trim(strg) 1719 read(strg1,*)(phfrq(mu,iph1l),mu=1,3*natom) 1720 else 1721 do nu=1,natom 1722 read(funit,*) (phfrq(((nu-1)*3)+mu,iph1l),mu=1,3) 1723 end do 1724 end if 1725 end if 1726 1727 if ((line(1:17)=='<dynamical_matrix')) then 1728 call rdfromline_value('dynamical_matrix',line,strg) 1729 if (strg/="") then 1730 ABI_MALLOC(work2,(3*natom,3*natom)) 1731 strg1=trim(strg) 1732 read(strg1,*) (work2(nu,1),nu=1,3*natom) 1733 do mu=2,3*natom-1 1734 read(funit,*)(work2(nu,mu),nu=1,3*natom) 1735 end do 1736 read(funit,'(a)',iostat=ios) readline 1737 call rmtabfromline(readline) 1738 line=adjustl(readline) 1739 call rdfromline_value('dynamical_matrix',line,strg) 1740 if (strg/="") then 1741 strg1=trim(strg) 1742 else 1743 strg1=trim(line) 1744 end if 1745 read(strg1,*) (work2(nu,3*natom),nu=1,3*natom) 1746 dynmat(1,:,:,:,:,iph1l) = reshape(work2,(/3,natom,3,natom/)) 1747 ABI_FREE(work2) 1748 else 1749 ABI_MALLOC(work2,(3*natom,3*natom)) 1750 do mu=1,3*natom 1751 read(funit,*)(work2(nu,mu),nu=1,3*natom) 1752 end do 1753 dynmat(1,:,:,:,:,iph1l) = reshape(work2,(/3,natom,3,natom/)) 1754 ABI_FREE(work2) 1755 end if 1756 end if 1757 1758 if ((line(1:8)==char(60)//char(47)//'phonon')) then 1759 iph1l = iph1l +1 1760 end if 1761 1762 if ((line(1:16)=='<strain_coupling')) then 1763 read(funit,'(a)',iostat=ios) readline 1764 call rdfromline("voigt",line,strg) 1765 strg1=trim(strg) 1766 read(strg1,*) voigt 1767 voigt = voigt + 1 ! 0 to 5 in the xml 1768 has_straincoupling = .true. 1769 irpt = 1 1770 end if 1771 1772 else 1773 ! Now treat the strain phonon coupling part 1774 if ((line(1:16)=='<strain_coupling')) then 1775 read(funit,'(a)',iostat=ios) readline 1776 call rdfromline("voigt",line,strg) 1777 strg1=trim(strg) 1778 read(strg1,*) voigt 1779 voigt = voigt + 1 ! 0 to 5 in the xml 1780 irpt = 1 1781 cycle 1782 end if 1783 1784 if(voigt>6)then 1785 write(message, '(4a)' )ch10,& 1786 & ' WARNING: the number of strain phonon coupling is superior to 6 in ',filename,ch10 1787 call wrtout(std_out,message,'COLL') 1788 exit 1789 end if 1790 1791 if ((line(1:22)=='<correction_force unit')) then 1792 call rdfromline_value('correction_force',line,strg) 1793 if (strg/="") then 1794 ABI_MALLOC(work2,(3,natom)) 1795 strg1=trim(strg) 1796 read(strg1,*) (work2(nu,1),nu=1,3) 1797 do mu=2,natom-1 1798 read(funit,*)(work2(nu,mu),nu=1,3) 1799 end do 1800 read(funit,'(a)',iostat=ios) readline 1801 call rmtabfromline(readline) 1802 line=adjustl(readline) 1803 call rdfromline_value('correction_force',line,strg) 1804 if (strg/="") then 1805 strg1=trim(strg) 1806 read(strg1,*) (work2(nu,natom),nu=1,3) 1807 else 1808 strg1=trim(line) 1809 read(strg1,*) (work2(nu,natom),nu=1,3) 1810 end if 1811 strain_coupling(voigt,:,:) = work2(:,:) 1812 ABI_FREE(work2) 1813 else 1814 ABI_MALLOC(work2,(3,natom)) 1815 do mu=1,natom 1816 read(funit,*)(work2(nu,mu),nu=1,3) 1817 end do 1818 strain_coupling(voigt,:,:) = work2(:,:) 1819 ABI_FREE(work2) 1820 end if 1821 end if 1822 1823 if ((line(1:11)=='<elastic3rd')) then 1824 call rdfromline_value('elastic3rd',line,strg) 1825 if (strg/="") then 1826 strg1=trim(strg) 1827 read(strg1,*) (elastic3rd(voigt,mu,1),mu=1,6) 1828 do nu=2,5 1829 read(funit,*) (elastic3rd(voigt,mu,nu),mu=1,6) 1830 end do 1831 else 1832 do nu=1,5 1833 read(funit,*) (elastic3rd(voigt,mu,nu),mu=1,6) 1834 end do 1835 end if 1836 read(funit,'(a)',iostat=ios) readline 1837 call rmtabfromline(readline) 1838 line=adjustl(readline) 1839 call rdfromline_value('elastic3rd',line,strg) 1840 if (strg/="") then 1841 strg1=trim(strg) 1842 read(strg1,*) (elastic3rd(voigt,mu,6),mu=1,6) 1843 else 1844 strg1=trim(line) 1845 read(strg1,*) (elastic3rd(voigt,mu,6),mu=1,6) 1846 end if 1847 has_anharmonics = .true. 1848 cycle 1849 end if 1850 1851 if ((line(1:29)=='<correction_strain_force unit')) then 1852 call rdfromline_value('correction_strain_force',line,strg) 1853 if (strg/="") then 1854 ABI_MALLOC(work2,(3*6,natom)) 1855 strg1=trim(strg) 1856 read(strg1,*) (work2(nu,1),nu=1,3*6) 1857 do mu=2,natom-1 1858 read(funit,*)(work2(nu,mu),nu=1,3*6) 1859 end do 1860 read(funit,'(a)',iostat=ios) readline 1861 call rmtabfromline(readline) 1862 line=adjustl(readline) 1863 call rdfromline_value('correction_strain_force',line,strg) 1864 if (strg/="") then 1865 strg1=trim(strg) 1866 read(strg1,*) (work2(nu,natom),nu=1,3*6) 1867 else 1868 strg1=trim(line) 1869 read(strg1,*) (work2(nu,natom),nu=1,3*6) 1870 end if 1871 elastic_displacement(voigt,:,:,:) = reshape(work2(:,:),(/6,3,natom/)) 1872 ABI_FREE(work2) 1873 else 1874 ABI_MALLOC(work2,(3*6,natom)) 1875 do mu=1,natom 1876 read(funit,*)(work2(nu,mu),nu=1,3*6) 1877 end do 1878 elastic_displacement(voigt,:,:,:) = reshape(work2(:,:),(/6,3,natom/)) 1879 ABI_FREE(work2) 1880 end if 1881 end if 1882 1883 if ((line(1:26)=='<correction_force_constant')) then 1884 found2=.false. 1885 do while (.not.found2) 1886 read(funit,'(a)',iostat=ios) readline 1887 call rmtabfromline(readline) 1888 line=adjustl(readline) 1889 if ((line(1:5)=='<data')) then 1890 call rdfromline_value('data',line,strg) 1891 if (strg/="") then 1892 ABI_MALLOC(work2,(3*natom,3*natom)) 1893 strg1=trim(strg) 1894 read(strg1,*) (work2(1,nu),nu=1,3*natom) 1895 do mu=2,3*natom-1 1896 read(funit,*)(work2(mu,nu),nu=1,3*natom) 1897 end do 1898 read(funit,'(a)',iostat=ios) readline 1899 call rmtabfromline(readline) 1900 line=adjustl(readline) 1901 call rdfromline_value('data',line,strg) 1902 if (strg/="") then 1903 strg1=trim(strg) 1904 read(strg1,*) (work2(3*natom,nu),nu=1,3*natom) 1905 else 1906 strg1=trim(line) 1907 read(strg1,*) (work2(3*natom,nu),nu=1,3*natom) 1908 end if 1909 phonon_strain(voigt)%atmfrc(:,:,:,:,irpt) = & 1910 & reshape(work2,(/3,natom,3,natom/)) 1911 ABI_FREE(work2) 1912 else 1913 ABI_MALLOC(work2,(3*natom,3*natom)) 1914 do mu=1,3*natom 1915 read(funit,*)(work2(mu,nu),nu=1,3*natom) 1916 end do 1917 phonon_strain(voigt)%atmfrc(:,:,:,:,irpt) =& 1918 & reshape(work2,(/3,natom,3,natom/)) 1919 ABI_FREE(work2) 1920 end if 1921 has_anharmonics = .true. 1922 end if 1923 if ((line(1:5)=='<cell')) then 1924 call rdfromline_value('cell',line,strg) 1925 if (strg/="") then 1926 strg1=trim(strg) 1927 read(strg1,*)(phonon_strain(voigt)%cell(mu,irpt),mu=1,3) 1928 else 1929 read(funit,*)(phonon_strain(voigt)%cell(mu,irpt),mu=1,3) 1930 end if 1931 irpt = irpt + 1 1932 found2=.true. 1933 cycle 1934 end if 1935 end do 1936 end if 1937 1938 if ((line(1:17)==char(60)//char(47)//'strain_coupling')) then 1939 ! set nrpt for the previous value of strain 1940 phonon_strain(voigt)%nrpt = irpt - 1 1941 ! restart the calculation of nrpt 1942 end if 1943 end if 1944 end if 1945 end do 1946 1947 1948 ! Reorder the ATMFRC 1949 ! Case 1: only local in the xml 1950 if (irpt1>0 .and. irpt2==0) then 1951 ifcs%cell(:,:) = int(cell_local(:,:)) 1952 ifcs%atmfrc(:,:,:,:,:) = local_atmfrc(:,:,:,:,:) 1953 ifcs%short_atmfrc(:,:,:,:,:) = local_atmfrc(:,:,:,:,:) 1954 ifcs%ewald_atmfrc(:,:,:,:,:) = zero 1955 1956 ! Case 2: only total in the xml 1957 else if(irpt1==0 .and. irpt2>0)then 1958 ifcs%cell(:,:) = int(cell_total(:,:)) 1959 ifcs%atmfrc(:,:,:,:,:) = total_atmfrc(:,:,:,:,:) 1960 ifcs%short_atmfrc(:,:,:,:,:) = zero 1961 ifcs%ewald_atmfrc(:,:,:,:,:) = total_atmfrc(:,:,:,:,:) 1962 1963 ! Case 3: local + total in the xml 1964 else if (irpt1>0 .and. irpt2>0)then 1965 if(irpt1 <= irpt2)then 1966 irpt3 = 0 1967 do ii=1,irpt2 1968 ifcs%cell(:,ii) = int(cell_total(:,ii)) 1969 ifcs%atmfrc(:,:,:,:,ii) = total_atmfrc(:,:,:,:,ii) 1970 do jj=1,irpt1 1971 if (all(abs(int(cell_local(:,jj))-ifcs%cell(:,ii))<tol16)) then 1972 ifcs%short_atmfrc(:,:,:,:,ii) = local_atmfrc(:,:,:,:,jj) 1973 irpt3 = irpt3 + 1 1974 end if 1975 end do 1976 end do 1977 if(irpt3 /= irpt1)then 1978 write(message, '(4a)' )ch10,& 1979 & ' There is several similar short IFC in ',filename,ch10 1980 ABI_BUG(message) 1981 end if 1982 else 1983 write(message, '(2a,I5,3a,I5,5a)' )ch10,& 1984 & ' The number of total IFC (',irpt2,') is inferior to ',ch10,& 1985 & ' the number of short range IFC (',irpt1,') in ',filename,ch10,& 1986 & ' This is not possible',ch10 1987 ABI_BUG(message) 1988 end if 1989 end if 1990 1991 ! Do some checks 1992 if (any(typat==0)) then 1993 write(message, '(a,a,a)' )& 1994 & ' Unable to read the type of atoms ',trim(filename),ch10 1995 ABI_ERROR(message) 1996 end if 1997 1998 if (any(abs(znucl)<tol16)) then 1999 write(message, '(a,a,a)' )& 2000 & ' Unable to read the atomic number ',trim(filename),ch10 2001 ABI_ERROR(message) 2002 end if 2003 2004 if (any(abs(all_amu)<tol16)) then 2005 write(message, '(a,a,a)' )& 2006 & ' Unable to read the atomic mass ',trim(filename),ch10 2007 ABI_ERROR(message) 2008 end if 2009 2010 close(unit=funit) 2011 2012 #endif 2013 2014 end if !End if master 2015 2016 !MPI BROADCAST 2017 call xmpi_bcast(energy,master, comm, ierr) 2018 call xmpi_bcast(all_amu,master, comm, ierr) 2019 call xmpi_bcast(dynmat,master, comm, ierr) 2020 call xmpi_bcast(elastic_constants,master, comm, ierr) 2021 call xmpi_bcast(epsilon_inf,master, comm, ierr) 2022 call xmpi_bcast(ifcs%nrpt,master, comm, ierr) 2023 call xmpi_bcast(ifcs%atmfrc,master, comm, ierr) 2024 call xmpi_bcast(ifcs%cell,master, comm, ierr) 2025 call xmpi_bcast(ifcs%ewald_atmfrc,master, comm, ierr) 2026 call xmpi_bcast(ifcs%short_atmfrc,master, comm, ierr) 2027 call xmpi_bcast(strain_coupling,master, comm, ierr) 2028 call xmpi_bcast(phfrq,master, comm, ierr) 2029 call xmpi_bcast(qph1l,master, comm, ierr) 2030 call xmpi_bcast(typat,master, comm, ierr) 2031 call xmpi_bcast(rprimd,master, comm, ierr) 2032 call xmpi_bcast(xcart,master, comm, ierr) 2033 call xmpi_bcast(zeff,master, comm, ierr) 2034 call xmpi_bcast(znucl,master, comm, ierr) 2035 do ii = 1,6 2036 call xmpi_bcast(phonon_strain(ii)%nrpt ,master, comm, ierr) 2037 call xmpi_bcast(phonon_strain(ii)%atmfrc ,master, comm, ierr) 2038 call xmpi_bcast(phonon_strain(ii)%cell ,master, comm, ierr) 2039 end do 2040 call xmpi_bcast(elastic3rd ,master, comm, ierr) 2041 call xmpi_bcast(elastic_displacement ,master, comm, ierr) 2042 call xmpi_bcast(has_anharmonics ,master, comm, ierr) 2043 2044 !Fill somes others variables 2045 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2046 call xcart2xred(natom,rprimd,xcart,xred) 2047 2048 !Re-generate symmetry operations from the lattice and atomic coordinates 2049 tolsym=tol8 2050 msym = 384 2051 ABI_MALLOC(spinat,(3,natom)) 2052 ABI_MALLOC(ptsymrel,(3,3,msym)) 2053 ABI_MALLOC(symafm,(msym)) 2054 ABI_MALLOC(symrel,(3,3,msym)) 2055 ABI_MALLOC(tnons,(3,msym)) 2056 use_inversion=1 2057 spinat = 0; 2058 symrel = 0; 2059 symafm = 0; 2060 tnons = 0 ; 2061 space_group = 0; 2062 call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym) 2063 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2064 call symfind(gprimd,msym,natom,nptsym,0,nsym,& 2065 & 0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred) 2066 2067 !Initialisation of crystal 2068 npsp = ntypat; timrev = 1 2069 ABI_MALLOC(title, (ntypat)) 2070 do ii=1,ntypat 2071 write(title(ii),'(a,i0)')"No title for typat ",ii 2072 end do 2073 2074 !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here 2075 call crystal_init(all_amu,Crystal,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,& 2076 & zion,znucl,timrev,.FALSE.,.FALSE.,title,& 2077 & symrel=symrel(:,:,1:nsym),tnons=tnons(:,1:nsym),symafm=symafm(1:nsym)) 2078 2079 !amu is not fill in crystal_init... 2080 Crystal%amu(:) = all_amu(:) 2081 2082 ABI_FREE(symrel) 2083 ABI_FREE(symafm) 2084 ABI_FREE(tnons) 2085 ABI_FREE(spinat) 2086 ABI_FREE(ptsymrel) 2087 2088 !if strcpling is set to 0 by the user, need to set the flag to false for 2089 !the initialisation of the effective potential 2090 if (present(strcpling))then 2091 if(strcpling == 0 )then 2092 has_anharmonics = .FALSE. 2093 end if 2094 end if 2095 2096 !Initialisation of eff_pot 2097 call effective_potential_init(crystal,eff_pot,energy,ifcs,ncoeff,nph1l,comm,& 2098 & dynmat=dynmat,elastic_constants=elastic_constants,& 2099 & elastic3rd=elastic3rd,elastic_displacement=elastic_displacement,& 2100 & epsilon_inf=epsilon_inf,strain_coupling=strain_coupling,& 2101 & phonon_strain=phonon_strain,phfrq=phfrq,qpoints=qph1l,& 2102 & has_anharmonicsTerms=has_anharmonics,zeff=zeff) 2103 2104 !DEALLOCATION OF ARRAYS 2105 ABI_FREE(all_amu) 2106 ABI_FREE(cell_local) 2107 ABI_FREE(cell_total) 2108 ABI_FREE(total_atmfrc) 2109 ABI_FREE(local_atmfrc) 2110 ABI_FREE(ifcs%atmfrc) 2111 ABI_FREE(ifcs%cell) 2112 ABI_FREE(ifcs%short_atmfrc) 2113 ABI_FREE(ifcs%ewald_atmfrc) 2114 ABI_FREE(dynmat) 2115 ABI_FREE(strain_coupling) 2116 ABI_FREE(phfrq) 2117 ABI_FREE(qph1l) 2118 ABI_FREE(title) 2119 ABI_FREE(typat) 2120 ABI_FREE(xcart) 2121 ABI_FREE(xred) 2122 ABI_FREE(zeff) 2123 ABI_FREE(zion) 2124 ABI_FREE(znucl) 2125 do ii = 1,6 2126 phonon_strain(ii)%nrpt = nrpt 2127 phonon_strain(ii)%atmfrc = zero 2128 phonon_strain(ii)%cell = 0 2129 ABI_FREE(phonon_strain(ii)%atmfrc) 2130 ABI_FREE(phonon_strain(ii)%cell) 2131 end do 2132 ABI_FREE(phonon_strain) 2133 ABI_FREE(elastic_displacement) 2134 2135 !DEALLOCATION OF TYPES 2136 call ifcs%free() 2137 call crystal%free() 2138 2139 end subroutine system_xml2effpot
m_effpot_xml/char_c2f [ Functions ]
NAME
char_c_to_f
FUNCTION
Helper function to convert a C string to a Fortran string Based on a routine by Joseph M. Krahn
INPUTS
c_string=C string
OUTPUT
f_string=Fortran string
SOURCE
3918 subroutine char_c2f(c_string,f_string) 3919 3920 use, intrinsic :: iso_c_binding, only : C_CHAR,C_NULL_CHAR 3921 !Arguments ------------------------------------ 3922 character(kind=C_CHAR,len=1),intent(in) :: c_string(*) 3923 character(len=*),intent(out) :: f_string 3924 !Local variables ------------------------------- 3925 integer :: ii 3926 !! ************************************************************************* 3927 ii=1 3928 do while(c_string(ii)/=C_NULL_CHAR.and.ii<=len(f_string)) 3929 f_string(ii:ii)=c_string(ii) ; ii=ii+1 3930 end do 3931 if (ii<len(f_string)) f_string(ii:)=' ' 3932 end subroutine char_c2f
m_effpot_xml/char_f2c [ Functions ]
NAME
char_f_to_c
FUNCTION
Helper function to convert a Fortran string to a C string Based on a routine by Joseph M. Krahn
INPUTS
f_string=Fortran string
OUTPUT
c_string=C string
SOURCE
3881 #if defined HAVE_XML 3882 3883 function char_f2c(f_string) result(c_string) 3884 3885 use, intrinsic :: iso_c_binding, only : C_CHAR,C_NULL_CHAR 3886 !Arguments ------------------------------------ 3887 character(len=*),intent(in) :: f_string 3888 character(kind=C_CHAR,len=1) :: c_string(len_trim(f_string)+1) 3889 !Local variables ------------------------------- 3890 integer :: ii,strlen 3891 !! ************************************************************************* 3892 strlen=len_trim(f_string) 3893 forall(ii=1:strlen) 3894 c_string(ii)=f_string(ii:ii) 3895 end forall 3896 c_string(strlen+1)=C_NULL_CHAR 3897 end function char_f2c