Project

General

Profile

cost_magfunctions.f90

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

Download (13 KB)

 
1
  Module cost_magfunctions 
2
      use CFML_GlobalDeps,                only: cp,sp,eps
3
      use CFML_crystallographic_symmetry, only: space_group_type
4

    
5
      use CFML_IO_Formats,                only: file_list_type
6
      use CFML_string_utilities,          only: l_case, setnum_std
7
      use CFML_Atom_TypeDef,              only: mAtom_List_Type
8
      use CFML_crystal_Metrics,           only: Crystal_Cell_Type
9
      use CFML_Keywords_Code_Parser,      only: NP_Max,NP_Refi,v_Vec,v_Shift,v_Bounds,v_BCon,v_Name,v_List, &
10
                                                VState_to_AtomsPar
11
      use CFML_Magnetic_Structure_Factors
12
      use CFML_Magnetic_Symmetry
13

    
14
      use prep_input
15

    
16
      implicit none
17
      private
18

    
19
      public::  General_cost_function, Readn_Set_CostFunctPars, Write_CostFunctPars, Write_FinalCost, &
20
                Write_SOL_mCFL
21

    
22
      logical,                      public :: err_cost=.false.
23
      character(len=132),           public :: err_mess_cost=" "
24
      type (space_group_type),      public :: SpG
25

    
26
      Integer, parameter,             public :: N_costf=1
27
      Integer,dimension(0:N_costf),   public :: Icost
28
      real,   dimension(0:N_costf),   public :: Wcost
29
      real,   dimension(0:N_costf),   public :: P_cost !Partial cost
30

    
31
      real,            public :: wavel
32
      character(len=3),public :: diff_mode="NUC"   ! XRA for x-rays, ELE for electrons
33

    
34
    Contains
35

    
36
!******************************************!
37
    Subroutine Readn_Set_CostFunctPars(file_cfl)
38
!******************************************!
39
       !---- Arguments ----!
40
       Type(file_list_type),   intent( in)  :: file_cfl
41
       !---- Local variables ----!
42
       character(len=132)   :: line
43
       real                 :: w,tol
44
       integer              :: i,ier,j,nrel
45

    
46
       Icost=0
47
       wcost=0.0
48

    
49
       i=0
50

    
51
       do j=1,file_cfl%nlines
52
          line=adjustl(file_cfl%line(j))
53
          line=l_case(line)
54
          if (len_trim(line) == 0) cycle
55
          if (line(1:1) =="!" .or. line(1:1) =="#") cycle
56
          i=index(line,"!")
57
          if( i /= 0) line=trim(line(1:i-1))
58

    
59
          select case (line(1:4))
60

    
61
              case ("opti")    !Optimization
62

    
63
              i=index(line,"f2mag")
64
              if(i == 0) then
65
                   Icost(0)=0; wcost(0)=0.0
66
                 else
67
                   Icost(0)=1
68
                   read(unit=line(i+9:),fmt=*,iostat=ier) w
69
                   if(ier /= 0) then
70
                      wcost(0)=1.0
71
                   else
72
                      wcost(0)= w
73
                   end if
74
                 end if
75

    
76
          end select
77

    
78
       end do
79

    
80
       return
81
    End Subroutine Readn_Set_CostFunctPars
82

    
83
!******************************************!
84
    Subroutine Write_CostFunctPars(lun)
85
!******************************************!
86
       !---- Arguments ----!
87
       integer,   intent( in)    :: lun
88
       !---- Local variables ----!
89
       character(len=132)   :: line
90
       real                 :: w,tol
91
       integer              :: i,ier
92

    
93

    
94
       Write(unit=lun,fmt="(/,a)")    "=================================="
95
       Write(unit=lun,fmt="(a)")      "= Cost functions to be minimized ="
96
       Write(unit=lun,fmt="(a,/)")    "=================================="
97

    
98
       do i=0,n_costf
99

    
100
          if(icost(i) == 0) cycle
101
          select case (i)
102

    
103
              case (0)    !Optimization "f2mag"
104
                 Write(unit=lun,fmt="(a,f8.4)") &
105
                 "  => Cost(F2mag): Optimization of Sum(|Obs-Scale*Calc|^2) / Sum (sigma^2)) / (Nobs-Nref), with weight: ",wcost(i)
106
          end select
107

    
108
       end do
109

    
110
       return
111
    End Subroutine Write_CostFunctPars
112

    
113
!******************************************!
114
    Subroutine Write_FinalCost(lun)
115
!******************************************!
116
       !---- Arguments ----!
117
       integer,   intent( in)    :: lun
118
       !---- Local variables ----!
119
       character(len=132)   :: line
120
       real                 :: w,tol
121
       integer              :: i
122

    
123

    
124
       Write(unit=lun,fmt="(/,a)")    "====================================="
125
       Write(unit=lun,fmt="(a)")      "= Final Cost of minimized functions ="
126
       Write(unit=lun,fmt="(a,/)")    "====================================="
127

    
128
       do i=0,n_costf
129

    
130
          if(icost(i) == 0) cycle
131

    
132
          select case (i)
133

    
134
              case (0)    !Optimization "F2mag"
135

    
136
                 Write(unit=lun,fmt="(a,f8.4,a,f12.4)") &
137
                 "  => Cost(F2mag): Optimization of Sum(|Obs-Scale*Calc|^2) / Sum(sigma^2)) / (Nobs-Nref), with weight: ",wcost(i),&
138
                 "  Final Cost: ",P_cost(0)
139

    
140
          end select
141

    
142
       end do
143

    
144
       return
145
    End Subroutine Write_FinalCost
146

    
147
!******************************************!
148
    Subroutine General_Cost_function(v,cost)
149
!******************************************!
150
      real,dimension(:),    intent( in):: v
151
      real,                 intent(out):: cost
152
      
153
       !---- Arguments ----!
154

    
155
      !---- Local variables ----!
156
      integer :: i,ic, nop, nlist=1, numv
157
      integer, dimension(1) :: List
158

    
159
      v_shift(1:NP_Refi)=v(1:NP_Refi)-v_vec(1:NP_Refi)  !Calculate the shifts w.r.t. old configuration
160
      v_vec(1:NP_Refi)=v(1:NP_Refi)
161

    
162
      numv=count(abs(v_shift) > 0.00001, dim=1)
163
      if(numv == 1) then
164
         i=maxloc(abs(v_shift),dim=1)
165
         List(1)=v_list(i)
166
         call VState_to_AtomsPar(mA,mode="S",MGp=MGp,Mag_dom=Mag_dom) !Update Atomic parameters with the proper constraints, 
167
         cost=0.0
168

    
169
         do ic=0,N_costf
170

    
171
            if(Icost(ic) == 0) cycle
172

    
173
            Select Case(ic)
174

    
175
               case(0)      !F2mag
176
                     call Calc_sqMiV_Data
177
                     call Cost_sqMiV(P_cost(0),Scale)
178
                     cost=cost+ P_cost(0)* WCost(0)
179

    
180
            End Select
181
         end do
182

    
183
      else            !New configuration
184

    
185
         call VState_to_AtomsPar(mA,mode="V",MGp=MGp,Mag_dom=Mag_dom)   !Update Atomic parameters with the proper constraints
186
         cost=0.0
187

    
188
         do ic=0,N_costf
189

    
190
            if(Icost(ic) == 0) cycle
191

    
192
            Select Case(ic)
193

    
194
               case(0)      !F2mag=sqMiV
195
                     call Calc_sqMiV_Data
196
                     call Cost_sqMiV(P_cost(0),Scale)
197
                     cost=cost+ P_cost(0)* WCost(0)
198

    
199
            End Select
200
         end do
201
      end if
202

    
203
      return
204
    End Subroutine General_Cost_function
205

    
206
!******************************************!
207
    Subroutine Cost_sqMiV(cost,Scale)
208
!******************************************!
209
       real,                 intent(out):: cost
210
       real,                 intent(in out):: Scale
211
       !---- Local variables ----!
212
       integer              :: j,n
213

    
214
       n=Oblist%Nobs
215
       Scale=sum( [(MhList%Mh(j)%sqMiV*Oblist%Ob(j)%Gobs*Oblist%Ob(j)%wGobs,j=1,n)])/ &
216
             sum( [(MhList%Mh(j)%sqMiV**2 * Oblist%Ob(j)%wGobs,j=1,n)] )
217
       cost=sum(([(Oblist%Ob(j)%wGobs* (Oblist%Ob(j)%Gobs-Scale*MhList%Mh(j)%sqMiV)**2,j=1,n)]))/(n-NP_Refi)
218
       return
219
    End Subroutine Cost_sqMiV
220

    
221
!******************************************!
222
Subroutine Write_SOL_mCFL(lun,file_cfl,mA,Mag_dom,comment)
223
!******************************************!
224
    !!----
225
    !!---- Subroutine Write_SOL_mCFL(lun,file_cfl,Cel,SpG,mA,Mag_dom,comment)
226
    !!----    integer,                  intent(in)      :: lun
227
    !!----    type(file_list_type),     intent (in out) :: file_cfl
228
    !!----    type (atom_list_type),    intent(in)      :: mA
229
    !!----    type (Magnetic_Domain_type),optional,intent(in)    :: Mag_dom
230
    !!----    character(len=*),optional,intent(in)      :: comment
231
    !!----
232
    !!----    Write a file mCFL
233
    !!----
234
    !!---- Created: January - 2012
235
    !!
236

    
237
       !---- Arguments ----!
238
       integer,                  intent(in)           :: lun
239
       type(file_list_type),     intent (in out)      :: file_cfl
240
       type (mAtom_list_type),   intent(in)           :: mA
241
       type (Magnetic_Domain_type),optional,intent(in):: Mag_dom
242
       character(len=*),optional,intent(in)           :: comment
243

    
244
       !----- Local variables -----!
245
       integer                         :: j,i,n,ier
246
       integer                         :: num_matom,num_skp,num_dom,ik,im,ip
247
       real,dimension(3)               :: Rsk,Isk
248
       real(kind=cp),dimension(12)     :: coef
249
       real(kind=cp)                   :: Ph
250
       character(len=132)              :: lowline,line
251
       character(len=30)               :: forma
252
       character(len=14)               :: pop
253
       logical                         :: skp_begin, bfcoef_begin, magdom_begin=.true.
254

    
255
       if(present(comment)) write(unit=lun,fmt="(a)") "TITLE "//trim(comment)
256
       write(unit=lun,fmt="(a)") "!  Automatically generated CFL file (Write_SOL_mCFL)"
257
       write(unit=lun,fmt="(a)") "!  "
258

    
259
      num_matom=0
260
      num_dom=0
261
      i=1
262

    
263
      do 
264
      i=i+1
265
      if(i >= file_cfl%nlines) exit
266

    
267
       lowline=l_case(adjustl(file_cfl%line(i)))
268
 
269
       if(lowline(1:6) == "magdom".and.magdom_begin) then
270
        num_dom=num_dom+1 
271
        ip=index(lowline,":")
272
        forma="(a  ,a14)"
273
        write(unit=forma(3:4),fmt="(i2)") ip
274
        write(pop,"(2f7.4)") Mag_Dom%Pop(1:2,num_dom)/sum(Mag_dom%Pop) !normalization to Sum_dom=1
275
        write(unit=file_cfl%line(i),fmt=forma) lowline(1:ip),pop
276

    
277
        do
278
         i=i+1
279
         lowline=adjustl(l_case(file_cfl%line(i)))         
280
         if(lowline(1:6) == "magdom") then
281
          num_dom=num_dom+1 
282
          ip=index(lowline,":")
283
          forma="(a  ,a14)"
284
          write(unit=forma(3:4),fmt="(i2)") ip
285
          write(pop,"(2f7.4)") Mag_Dom%Pop(1:2,num_dom)/sum(Mag_dom%Pop) !normalization to Sum_dom=1
286
          write(unit=file_cfl%line(i),fmt=forma) lowline(1:ip),pop
287
         else
288
          i=i-1
289
          magdom_begin=.false.
290
          exit
291
         endif
292
        enddo
293
       endif! end magdom
294

    
295
       if(lowline(1:5) == "matom") then
296
          num_matom=num_matom+1 !max NMatom=mA%natoms
297
          num_skp=0
298
          skp_begin=.true.
299
          bfcoef_begin=.true.
300
          write(unit=lun,fmt="(a)") trim(file_cfl%line(i))
301
          cycle
302
       endif
303

    
304
       if(lowline(1:3) == "skp".and.skp_begin) then
305
        num_skp=num_skp+1 !max mA%atom(num_matom)%nvk 
306
        read(unit=lowline(4:),fmt=*,iostat=ier) ik,im,Rsk,Isk,Ph
307
        if(MGp%Sk_type == "Spherical_Frame") then
308
         Rsk(:)=mA%atom(num_matom)%Spher_Skr(:,ik)
309
         Isk(:)=mA%atom(num_matom)%Spher_Ski(:,ik)
310
         Ph=mA%atom(num_matom)%mphas(ik)
311
        else
312
         Rsk(:)=mA%atom(num_matom)%Skr(:,ik)
313
         Isk(:)=mA%atom(num_matom)%Ski(:,ik)
314
         Ph=mA%atom(num_matom)%mphas(ik)        
315
        endif
316
        
317
        write(unit=file_cfl%line(i),fmt='(a,i8,i3,7f8.3)') 'skp',ik,im,Rsk,Isk,Ph
318
        do
319
         i=i+1
320
         lowline=adjustl(l_case(file_cfl%line(i)))         
321
         if(lowline(1:3) == "skp") then
322
          num_skp=num_skp+1
323
          read(unit=lowline(4:),fmt=*,iostat=ier) ik,im,Rsk,Isk,Ph
324
          if(MGp%Sk_type == "Spherical_Frame") then
325
           Rsk(:)=mA%atom(num_matom)%Spher_Skr(:,ik)
326
           Isk(:)=mA%atom(num_matom)%Spher_Ski(:,ik)
327
           Ph=mA%atom(num_matom)%mphas(ik)
328
          else
329
           Rsk(:)=mA%atom(num_matom)%Skr(:,ik)
330
           Isk(:)=mA%atom(num_matom)%Ski(:,ik)
331
           Ph=mA%atom(num_matom)%mphas(ik)        
332
          endif
333
          write(unit=file_cfl%line(i),fmt='(a,i8,i3,7f8.3)') 'skp',ik,im,Rsk,Isk,Ph
334
         else
335
          i=i-1
336
          skp_begin=.false.
337
          exit
338
         endif
339
        enddo
340

    
341
       endif! end Rsk,Isk,Ph
342

    
343
       forma="(a6,i8,i3,  f8.3)"
344

    
345
       if(lowline(1:6) == "bfcoef".and.bfcoef_begin) then
346
        num_skp=num_skp+1 !max mA%atom(num_matom)%nvk 
347
        read(unit=lowline(7:),fmt=*,iostat=ier) ik,im
348
        n=abs(MGp%nbas(im))
349
        write(unit=forma(11:12),fmt="(i2)") n+1
350
        read(unit=lowline(7:),fmt=*,iostat=ier) ik,im,coef(1:n),ph
351
        coef(1:n)= mA%atom(num_matom)%cbas(1:n,ik) 
352
        Ph=mA%atom(num_matom)%mphas(ik)        
353
        write(unit=file_cfl%line(i),fmt=forma) 'bfcoef',ik,im,coef(1:n),Ph        
354
        do
355
         i=i+1
356
         lowline=adjustl(l_case(file_cfl%line(i)))         
357
         if(lowline(1:6) == "bfcoef") then
358
          num_skp=num_skp+1
359
          read(unit=lowline(7:),fmt=*,iostat=ier) ik,im
360
          n=abs(MGp%nbas(im))
361
          write(unit=forma(11:12),fmt="(i2)") n+1
362
          read(unit=lowline(7:),fmt=*,iostat=ier) ik,im,coef(1:n),ph
363
          coef(1:n)= mA%atom(num_matom)%cbas(1:n,ik) 
364
          Ph=mA%atom(num_matom)%mphas(ik)        
365
          write(unit=file_cfl%line(i),fmt=forma) 'bfcoef',ik,im,coef(1:n),Ph
366
         else
367
          i=i-1
368
          bfcoef_begin=.false.
369
          exit
370
        endif
371
        enddo
372
       endif! end bfcoef
373

    
374
       write(unit=lun,fmt="(a)") trim(file_cfl%line(i))
375

    
376
      end do
377

    
378
    End Subroutine Write_SOL_mCFL
379

    
380
  End Module cost_magfunctions