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 xrays, 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:i1))

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(ObsScale*Calc^2) / Sum (sigma^2)) / (NobsNref), 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(ObsScale*Calc^2) / Sum(sigma^2)) / (NobsNref), 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)%GobsScale*MhList%Mh(j)%sqMiV)**2,j=1,n)]))/(nNP_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=i1

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=i1

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=i1

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
