Project

General

Profile

prep_input.f90

Zaharko, Oksana, 05/02/2012 12:40 PM

Download (12.6 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
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, Copy_MagDom_List, &
18
           Write_ObsCalc, Mag_Dom, mA
19

    
20
 integer, save               :: lun=1,lan=2,lin=3
21
 integer, save               :: Nf2
22
 real, public                :: Scale
23

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

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

    
65
 type (file_list_type)             :: file_cfl
66
 type (MagSymm_k_Type)             :: MGp
67
 type (MagHD_Type)                 :: Mh
68
 type (MagHD_List_Type)            :: Mhlist
69
 type (mObservation_Type)          :: Ob
70
 type (mObservation_List_Type)     :: Oblist
71
 type (Crystal_Cell_Type)          :: Cell
72
 type (mAtom_List_Type),save       :: mA,mA_clone
73
 type (Magnetic_Domain_type),save  :: Mag_Dom,Mag_dom_clone
74
 
75
contains
76

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

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

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

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

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

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

    
128
1    rewind(unit=lin)
129

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

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

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

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

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

    
184
2    close(lin)     
185

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

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

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

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

    
217
!******************************************!
218
    Subroutine Calc_sqMiV_Data(lun)
219
!******************************************!
220
! gets magn str Factor and Miv from Calc_Magnetic_StrF_MiV
221
! 
222
       !---- Argument ----!
223
       integer, optional,           intent(in) :: lun
224

    
225
       !---- Local variables ----!
226
       integer                     :: j,n
227
       real                        :: cost
228
       character(len=1)            :: mode
229

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

    
234
         !Calculate magnetic structure factor and magnetic interaction vector
235
         ! as mode='y' MiV w.r.t. cartesian crystallographic frame
236
         mode='y'
237
         Mh=MhList%Mh(j)
238
!         call Calc_Magnetic_StrF_MiV(Cell,MGp,mA,Mh) ! was withput domains 
239
         call Calc_Magnetic_StrF_MiV_Dom(Cell,MGp,mA,Mag_Dom,Mh) !CFML_Msfac
240

    
241
!----    without mode components with respect to direct cell system {e1,e2,e3}
242
!----    with mode components with respect to the cartesian frame
243
         MhList%Mh(j)%sqMiV=Mh%sqMiV
244
       end do !end loop over reflections
245

    
246
!---     Cost_sqMiV(cost,Scale)
247
       if (present(lun)) then
248
          Scale=sum( [(MhList%Mh(j)%sqMiV*Oblist%Ob(j)%Gobs *Oblist%Ob(j)%wGobs,j=1,n)])/ &
249
                sum( [(MhList%Mh(j)%sqMiV**2 *Oblist%Ob(j)%wGobs,j=1,n)] )
250

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

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

    
261
       end if
262

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

    
284
       n=mA%natoms
285
       call Allocate_mAtom_List(n,mA_clone)
286
       mA_clone%atom(1:n)=mA%atom(1:n)
287
       return
288
    End Subroutine Copy_mAtom_List
289

    
290
!******************************************!
291
    !!----
292
    !!---- Subroutine Copy_MagDom_List(Mag_dom, Mag_dom_clone)
293
    !!----    type(mAtom_list_type),   intent(in)  :: Mag_dom      !  In -> Atom list
294
    !!----    type(mAtom_list_type),   intent(out) :: Mag_dom_clone!   Out -> Atom list
295
    !!----
296
    !!----    Subroutine to copy an mag domain list to another one
297
    !!----
298
    !!---- Created: January - 2012
299
    !!
300
!******************************************!
301
    Subroutine Copy_MagDom_List(Mag_dom, Mag_dom_clone)
302
!******************************************!
303
       !---- Arguments ----!
304
       type(Magnetic_Domain_type),   intent(in)  :: Mag_dom  
305
       type(Magnetic_Domain_type),   intent(out) :: Mag_dom_clone
306
       !---- Local Variables ----!
307

    
308
       
309
       Mag_dom_clone%nd=Mag_dom%nd
310
       Mag_dom_clone%Chir=Mag_dom%Chir
311
       Mag_dom_clone%DMat=Mag_dom%DMat
312
       Mag_dom_clone%Pop=Mag_dom%Pop
313
       Mag_dom_clone%Mpop=Mag_dom%Mpop
314
       Mag_dom_clone%Lpop=Mag_dom%Lpop
315
       Mag_dom_clone%Lab=Mag_dom%Lab
316

    
317
       return
318
    End Subroutine Copy_MagDom_List
319

    
320
!******************************************!
321
    Subroutine Write_ObsCalc(file_cfl)
322
!******************************************!
323

    
324
       !---- Arguments ----!
325
       type(file_list_type),                intent (in)     :: file_cfl
326
       !---- Local variables ----!
327
       character(len=132)   :: line
328
       integer              :: i,j
329

    
330
        i=0
331
        do j=1,file_cfl%nlines
332
          line=adjustl(file_cfl%line(j))
333
          line=l_case(line)
334
          if (line(1:1) ==" ") cycle
335
          if (line(1:1) =="!") cycle
336
          i=index(line,"!")
337

    
338
          if( i /= 0) line=trim(line(1:i-1))
339

    
340
          select case (line(1:4))
341
          case ("opti")    !Optimization
342

    
343
          i=index(line,"f2mag")
344
           if(i /= 0) then
345
            write(lun,'(a)') '     hkl                      Iobs      Ical'
346
            do i=1,Nf2
347
             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
348
            enddo
349
           endif
350

    
351
          end select
352
        enddo
353
   
354
   End Subroutine Write_ObsCalc
355
   
356
end module prep_input