Project

General

Profile

prep_input.f90

Zaharko, Oksana, 22/01/2012 10:29 AM

Download (11.3 KB)

 
1
Module prep_input
2
 !---- Use Modules ----! 
3
 use CFML_GlobalDeps,             only: sp, cp, eps
4
 use CFML_Atom_TypeDef,           only: mAtom_list_Type, allocate_mAtom_list
5
 use CFML_String_Utilities,       only: l_case
6
 use CFML_Reflections_Utilities,  only: Hkl_s
7
 use CFML_IO_Formats,             only: file_list_type
8
 use CFML_Propagation_vectors,    only: K_Equiv_Minus_K
9
 use CFML_crystal_Metrics,        only: Crystal_Cell_Type
10
 use CFML_Keywords_Code_Parser,   only: NP_Refi
11
 use CFML_Magnetic_Symmetry 
12
 use CFML_Magnetic_Structure_Factors
13

    
14
 implicit none
15

    
16
   !---- List of public subroutines ----!
17
 public :: Readn_Set_sqMiV_Data, Readn_Set_Job, Calc_sqMiV_Data, Copy_mAtom_List, Write_ObsCalc
18

    
19
 integer, save               :: lun=1,lan=2,lin=3
20
 integer, save               :: Nf2
21
 real, public                :: Scale
22
 
23
   !---- Other types are given belove the new types
24
   
25
    !!----
26
    !!---- TYPE :: mOBSERVATION_TYPE
27
    !!--..
28
    !!---- Type, public :: mObservation_Type
29
    !!----    integer                 :: ncont ! number of contributing reflections
30
    !!-----   integer,dimension(3,12) :: Hd    ! Indices of contributing reflections
31
    !!-----   integer,dimension(  12) :: icod  ! Pointer to the number of domain
32
    !!-----   integer,dimension(  12) :: p     ! Pointer to the list of independent reflections
33
    !!----    real(kind=sp)           :: Gobs  ! Observed  Magnetic Intensity
34
    !!----    real(kind=sp)           :: SGobs ! Sigma of  Gobs
35
    !!---- End Type mObservation_Type
36
    !!----
37
    !!---- Created: January - 2012
38
    !!
39
    Type, public :: mObservation_Type
40
       integer                                  :: ncont ! number of contributing reflections
41
       real(kind=sp),allocatable,dimension(:,:) :: Hd    ! (3,ncont) Indices of contributing reflections
42
       integer,allocatable,dimension(:)         :: icod  ! (ncont)   Pointer to the number of domain
43
       integer,allocatable,dimension(:)         :: p     ! (ncont)   Pointer to the list of independent reflections
44
       real(kind=sp)                            :: Gobs  ! Observed  Magnetic Intensity
45
       real(kind=sp)                            :: sGobs ! Sigma of  Gobs
46
       real(kind=sp)                            :: wGobs ! Weight of Gobs, usually 1/sGobs^2
47
    End Type mObservation_Type
48

    
49
    !!----
50
    !!---- TYPE :: mOBSERVATION_LIST_TYPE
51
    !!--..
52
    !!---- Type, public :: mObservation_List_Type
53
    !!----    integer                                          :: Nobs  ! Number of Observations
54
    !!----    type(mobservation_type),allocatable,dimension(:) :: Ob    ! Observation (F2, etc...)
55
    !!---- End Type mObservation_List_Type
56
    !!----
57
    !!---- Created: January - 2012
58
    !!
59
    Type, public :: mObservation_List_Type
60
       integer                                          :: Nobs  ! Number of Observations
61
       type(mObservation_type),allocatable,dimension(:) :: Ob    ! Observation (F2, etc...)
62
    End Type mObservation_List_Type
63

    
64
 type (file_list_type)             :: file_cfl
65
 type (MagSymm_k_Type)             :: MGp
66
 type (MagH_Type)                  :: Mh
67
 type (MagH_List_Type)             :: Mhlist
68
 type (mAtom_list_type)            :: mA,mA_clone
69
 type (mObservation_Type)          :: Ob
70
 type (mObservation_List_Type)     :: Oblist
71
 type (Crystal_Cell_Type)          :: Cell
72
 
73
contains
74

    
75
!---- Subroutines ----!
76
!******************************************!
77
    Subroutine Readn_Set_sqMiV_Data(file_cfl)
78
!******************************************!
79
!    reads integrated intensities for the list of reflections given in *.int
80
       !---- Arguments ----!
81
       type(file_list_type),                intent (in)     :: file_cfl
82

    
83
       !---- Local variables ----!
84
       character(len=256)          :: title,formt
85
       character(len=256)          :: line,fildat
86
       character(len=1)            :: sig
87
       integer                     :: ier,num_k,m,ih,ik,il,ncont,iv,i,j,ityp,ipow
88
       real,   dimension(3)        :: vk,v_k
89
       real                        :: Gobs,SGobs,lambda
90
       logical                     :: esta
91
    
92
          esta=.false.
93
          do i=1,file_cfl%nlines
94
            line=adjustl(l_case(file_cfl%line(i)))
95
            if(line(1:10) == "f2mag_data") then
96
              line=file_cfl%line(i)
97
              j=index(line,"!")
98
              if(j /= 0) line=trim(line(1:j-1))
99
              fildat= adjustl(line(12:))
100
              inquire(file=trim(fildat),exist=esta)
101
              exit
102
            end if
103
          end do
104

    
105
          if(.not. esta) then
106
            write(unit=*,fmt="(a)") " => No .int has been provided! "
107
          end if
108

    
109
     open(unit=lin,file=trim(fildat), status="old",action="read")
110
! temporarily programmed for the special case k,-k, i.e. m=1,2     
111
     read(unit=lin,fmt=*) title
112
     read(unit=lin,fmt=*) formt
113
     read(unit=lin,fmt=*) lambda,ityp,ipow
114
     read(unit=lin,fmt=*) num_k
115

    
116
     do i=1,num_k
117
      read(unit=lin,fmt=*) m,vk
118
     enddo
119

    
120
     Nf2=0
121
     do !just to get Nf2
122
      read(unit=lin,fmt=*,end=1) ih,ik,il,m,Gobs,SGobs,ncont
123
      Nf2=Nf2+1
124
     enddo
125

    
126
1    rewind(unit=lin)
127

    
128
     if(allocated(Oblist%Ob)) deallocate(Oblist%Ob)
129
     allocate(Oblist%Ob(Nf2))
130
     
131
     if(allocated(Mhlist%Mh)) deallocate(Mhlist%Mh)
132
     allocate(Mhlist%Mh(Nf2))
133
     Mhlist%Nref= Nf2
134

    
135
! temporarily programmed for twin coefficient ncont=1
136
     do i=1,Nf2
137
      allocate(Oblist%Ob(i)%Hd(1:3,1))
138
     enddo
139
     Oblist%Nobs= Nf2
140
     
141
     read(unit=lin,fmt=*) title
142
     read(unit=lin,fmt=*) formt
143
     read(unit=lin,fmt=*) lambda,ityp,ipow
144
     read(unit=lin,fmt=*) num_k
145

    
146
     do i=1,num_k
147
      if(i.eq.1) then
148
       read(unit=lin,fmt=*) m,vk
149
      else
150
       read(unit=lin,fmt=*) m,v_k
151
       if(maxval(abs(vk+v_k)).gt.eps) then
152
        write(unit=*,fmt="(a)") " => Only k and -k are handled presently! "
153
        stop
154
       endif
155
      endif
156
     enddo
157

    
158
     do i=1,Nf2
159
      read(unit=lin,fmt=*,end=2) ih,ik,il,m,Oblist%Ob(i)%Gobs,Oblist%Ob(i)%SGobs,Oblist%Ob(i)%ncont
160
      if(m==1) then
161
       Oblist%Ob(i)%Hd(:,1)=real((/ih,ik,il/)) + vk(:)
162
       j=1
163
      endif
164
      if(m==2) then
165
       Oblist%Ob(i)%Hd(:,1)=real((/ih,ik,il/)) + v_k(:)
166
       j=-1
167
      endif
168

    
169
      !construct partially the object Mh
170
         if(j== 1) sig='+'
171
         if(j==-1) sig='-'
172
         Mh%signp=real(-j)  ! sign "+" for H-k and "-" for H+k
173
         iv=abs(j)
174
         Mh%num_k=iv
175
         Mh%h= real((/ih,ik,il/)) - Mh%signp*MGp%kvec(:,iv)
176
         Mh%s = hkl_s(Mh%h,Cell)
177
         Mh%keqv_minus=K_Equiv_Minus_K(MGp%kvec(:,iv),MGp%latt)
178
         MhList%Mh(i)=Mh
179
     enddo
180
       Oblist%Ob(:)%wGobs=1.d0/max([(Oblist%Ob(i)%sGobs**2,i=1,Nf2)],eps)
181

    
182
2    close(lin)     
183

    
184
    End Subroutine Readn_Set_sqMiV_Data
185
    
186
!******************************************!
187
    Subroutine Readn_Set_Job(file_cfl)
188
!******************************************!
189
       !---- Arguments ----!
190
       Type(file_list_type),   intent( in)  :: file_cfl
191

    
192
       !---- Local variables ----!
193
       character(len=132)   :: line
194
       integer              :: i,j
195

    
196
       i=0
197
       do j=1,file_cfl%nlines
198
          line=adjustl(file_cfl%line(j))
199
          line=l_case(line)
200
          if (line(1:1) ==" ") cycle
201
          if (line(1:1) =="!") cycle
202
          i=index(line,"!")
203
          if( i /= 0) line=trim(line(1:i-1))
204

    
205
          select case (line(1:4))
206
          case ("opti")    !Optimization
207
              i=index(line,"f2mag")
208
              if(i /= 0) then
209
               call Calc_sqMiV_Data(lun)
210
              endif
211
          end select
212
       end do
213
    End Subroutine Readn_Set_Job 
214

    
215
!******************************************!
216
    Subroutine Calc_sqMiV_Data(lun)
217
!******************************************!
218
! gets magn str Factor and Miv from Calc_Magnetic_StrF_MiV
219
! 
220
       !---- Argument ----!
221
       integer, optional,           intent(in) :: lun
222
       !---- Local variables ----!
223
       integer                     :: j,n
224
       real                        :: cost
225
       character(len=1)            :: mode
226

    
227
       if (present(lun)) write(lun,'(a)') '   hkl                      Imag   sImag   Icalc'
228
       n=Oblist%Nobs       
229
       do j=1,Nf2 ! This is loop over reflections with Int
230

    
231
         !Calculate magnetic structure factor and magnetic interaction vector
232
         ! as mode='y' MiV w.r.t. cartesian crystallographic frame
233
         mode='y'
234
         Mh=MhList%Mh(j)
235
         call Calc_Magnetic_StrF_MiV(Cell,MGp,mA,Mh)
236
!----    without mode components with respect to direct cell system {e1,e2,e3}
237
!----    with mode components with respect to the cartesian frame
238
         MhList%Mh(j)%sqMiV=Mh%sqMiV
239
       end do !end loop over reflections
240

    
241
!---     Cost_sqMiV(cost,Scale)
242
       if (present(lun)) then
243
          Scale=sum( [(MhList%Mh(j)%sqMiV*Oblist%Ob(j)%Gobs *Oblist%Ob(j)%wGobs,j=1,n)])/ &
244
                sum( [(MhList%Mh(j)%sqMiV**2 *Oblist%Ob(j)%wGobs,j=1,n)] )
245

    
246
         do j=1,Nf2 ! This is loop over reflections with Int
247
          write(unit=lun,fmt="(6f8.4)") Oblist%Ob(j)%Hd(:,1),&
248
                                        Oblist%Ob(j)%Gobs,Oblist%Ob(j)%SGobs,Scale*MhList%Mh(j)%sqMiV
249
         end do !end loop over reflections
250
         cost=sum( ([(Oblist%Ob(j)%wGobs* (Oblist%Ob(j)%Gobs-Scale*MhList%Mh(j)%sqMiV)**2,j=1,n)]) )/(n-NP_Refi)
251

    
252
         write(unit=lun,fmt="(/,a)")    "=================================="
253
         write(unit=lun,fmt="(a,f12.0)") 'Initial Cost(F2mag): Sum(|Obs-Scale*Calc|^2) / Sum(sigma^2) ) / (Nobs-Nref) =',cost
254
         write(unit=lun,fmt="(a,/)")    "=================================="
255

    
256
       end if
257

    
258
    End Subroutine Calc_sqMiV_Data
259
   
260
!******************************************!
261
    !!----
262
    !!---- Subroutine Copy_mAtom_List(mA, mA_clone)
263
    !!----    type(mAtom_list_type),   intent(in)  :: mA      !  In -> Atom list
264
    !!----    type(mAtom_list_type),   intent(out) :: mA_clone     !   Out -> Atom list
265
    !!----
266
    !!----
267
    !!----    Subroutine to copy an atom list to another one
268
    !!----
269
    !!---- Created: January - 2012
270
    !!
271
!******************************************!
272
    Subroutine Copy_mAtom_List(mA, mA_clone)
273
!******************************************!
274
       !---- Arguments ----!
275
       type(mAtom_list_type),   intent(in)  :: mA  
276
       type(mAtom_list_type),   intent(out) :: mA_clone
277
       !---- Local variables ----!
278
       integer                    :: n
279

    
280
       n=mA%natoms
281
       call Allocate_mAtom_List(n,mA_clone)
282
       mA_clone%atom(1:n)=mA%atom(1:n)
283
       return
284
    End Subroutine Copy_mAtom_List
285

    
286
!******************************************!
287
    Subroutine Write_ObsCalc(file_cfl)
288
!******************************************!
289

    
290
       !---- Arguments ----!
291
       type(file_list_type),                intent (in)     :: file_cfl
292
       !---- Local variables ----!
293
       character(len=132)   :: line
294
       integer              :: i,j
295

    
296
        i=0
297
        do j=1,file_cfl%nlines
298
          line=adjustl(file_cfl%line(j))
299
          line=l_case(line)
300
          if (line(1:1) ==" ") cycle
301
          if (line(1:1) =="!") cycle
302
          i=index(line,"!")
303

    
304
          if( i /= 0) line=trim(line(1:i-1))
305

    
306
          select case (line(1:4))
307
          case ("opti")    !Optimization
308

    
309
          i=index(line,"f2mag")
310
           if(i /= 0) then
311
            write(lun,'(a)') '     hkl                      Iobs      Ical'
312
            do i=1,Nf2
313
             write(unit=lun,fmt='(3(f8.4),2x,2(f12.2))') Oblist%Ob(i)%Hd(:,1),Oblist%Ob(i)%Gobs,Scale*Mhlist%Mh(i)%sqMiV
314
            enddo
315
           endif
316

    
317
          end select
318
        enddo
319
   
320
   End Subroutine Write_ObsCalc
321
   
322
end module prep_input