1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
|
!
! Copyright (C) 2017 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------
PROGRAM ibrav2cell
!----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : pi, ANGSTROM_AU
!USE powell, ONLY : POWELL_MIN
USE lmdif_module, ONLY : lmdif1
!
IMPLICIT NONE
CHARACTER(len=1024) :: line
INTEGER,PARAMETER :: npar = 6+3 ! celldm(6) and 3 angles
INTEGER,PARAMETER :: nibrav = 20
INTEGER,PARAMETER :: ibrav_list(nibrav) = (/1,2,3,-3,4,5,-5,6,7,8,9,-9,91,10,11,12,-12,13,-13,14/)
INTEGER :: ibrav, ios, ii, i, info
REAL(DP) :: celldm(6), angle(3), alat, chisq, chisq_aux
!
REAL(DP) :: at(3,3), omega, R(3,3), par(npar), par_aux(npar), vchisq(npar)
REAL(DP),PARAMETER :: grad_to_rad = pi/180
! REAL(DP) :: xi(npar,npar)
INTEGER :: lwa, iwa(npar)
REAL(DP),ALLOCATABLE :: wa(:)
!
LOGICAL,EXTERNAL :: matches
WRITE(*,*) "Enter the unit of measur (angstrom, bohr) or alat in bohr units, or alat in Angstrom units, followed by ' A'"
READ(*,"(a1024)") line
IF(matches("angstrom",line)) THEN
alat=ANGSTROM_AU
ELSE IF(matches("bohr",line)) THEN
alat=1._dp
ELSE
READ(line,*,iostat=ios) alat
IF(ios/=0) CALL errore("scan_ibrav","Could not understand alat '"//TRIM(line)//"'",1)
IF(matches("A",line)) alat=alat*ANGSTROM_AU
ENDIF
WRITE(*,'("alat (bohr)",f12.6)') alat
WRITE(*,*) "Enter the cell basis vectors (one per line)"
DO i = 1,3
READ(*,*) at(:,i)
ENDDO
at = at*alat
WRITE(*,'("Requested axes in Bohr units:")')
WRITE(*, '("at1", 6f14.6)') at(:,1)
WRITE(*, '("at2", 6f14.6)') at(:,2)
WRITE(*, '("at3", 6f14.6)') at(:,3)
lwa = npar**2 + 6*npar
ALLOCATE(wa(lwa))
DO ii = 1, nibrav
ibrav = ibrav_list(ii)
!xi = 0._dp
!FORALL(i=1:npar) xi(i,i) = 1._dp
par(1) = SQRT(SUM(at**2))/3
par(2) = 1._dp
par(3) = 1._dp
par(4) = 0.1_dp
par(5) = 0.1_dp
par(6) = 0.1_dp
par(7) = 0._dp
par(8) = 0._dp
par(9) = 0._dp
WRITE(*,'("Scanning ibrav ",i3,"")') ibrav
CALL lmdif1(optimize_this_s, npar, npar, par, vchisq, 1.d-15, info, iwa, wa, lwa)
IF(info>0 .and. info<5) THEN
PRINT*, "Minimization succeeded"
ELSEIF(info>=5) THEN
PRINT*, "Minimization stopped before convergence"
ELSEIF(info<=0) THEN
PRINT*, "Minimization error", info
STOP
ENDIF
chisq = vchisq(1)
IF(chisq<1.d-6)THEN
WRITE(*,'(/,/,"____________ MATCH (chisq=",g7.1,") ____________")') chisq
WRITE(*, '(" ibrav = ",i3)') ibrav
DO i = 1,6
par_aux = par
par_aux(i) = par_aux(i) * .5_dp
!chisq_aux = optimize_this(par_aux)
CALL optimize_this_s(npar, npar, par_aux, vchisq, info)
chisq_aux = vchisq(1)
IF(chisq_aux/=chisq)THEN
WRITE(*, '(" celldm(",i2,") = ", f14.9)') i,par(i)
ENDIF
ENDDO
IF(ANY(ABS(par(7:9))>1.d-12))THEN
WRITE(*,'(a,/,a)') "WARNING! Cell is rotated. Atomic positions will also need",&
" to be rotated if they are not in crystal coords!"
WRITE(*, '("angles (around x,y,z)", 6f14.3)') par(7:9)
ENDIF
WRITE(*, '("at1", 6f14.6)') at(:,1)
WRITE(*, '("at2", 6f14.6)') at(:,2)
WRITE(*, '("at3", 6f14.6)') at(:,3)
WRITE(*,*)
ENDIF
ENDDO
!
CONTAINS
SUBROUTINE optimize_this_s(m_,n_,p_,f_,i_)
IMPLICIT NONE
INTEGER,INTENT(in) :: m_, n_
INTEGER,INTENT(inout) :: i_
REAL(DP),INTENT(out) :: f_(m_)
REAL(DP),INTENT(in) :: p_(n_)
REAL(DP) :: celldm_(6), angle_(3), at_(3,3), R(3,3), omega_, pars_(n_), penality
INTEGER :: ierr
CHARACTER(len=32) :: errormsg
! Global variables from main function: ibrav, at
pars_ = p_ ! p_ is read only!
!WRITE(*, '("A:", 6f10.4,3x,3f10.4)') pars_
CALL check_bounds(pars_, penality)
!WRITE(*, '("B:", 6f10.4,3x,3f10.4)') pars_
celldm_ = pars_(1:6)
angle_ = pars_(7:9)*grad_to_rad
!WRITE(*, '(6(f14.4,2x),3x,3f10.4)') celldm_, angle_
CALL latgen_lib( ibrav, celldm_, at_(:,1), at_(:,2), at_(:,3), omega_, ierr, errormsg)
!IF(ierr/=0) penalty=penalty*10
!
IF (ANY(angle_/=0._dp)) THEN
R = rot(angle_(1), angle_(2), angle_(3))
at_ = matmul(R,at_)
ENDIF
f_ = 0._dp
f_(1) = SUM( (at-at_)**2 )*penality
END SUBROUTINE
SUBROUTINE check_bounds(pars_, penalty)
IMPLICIT NONE
REAL(DP),INTENT(inout) :: pars_(npar), penalty
REAL(DP),PARAMETER :: infty = 1.d+100, eps=1.d-6
REAL(DP),PARAMETER :: par_min(npar) = (/ eps, eps, eps, -.5_dp+eps, -.5_dp+eps, -1._dp, -180._dp, -180._dp, -180._dp /)
REAL(DP),PARAMETER :: par_max(npar) = (/ infty, infty, infty, 1._dp-eps, 1._dp-eps, 1._dp-eps, 180._dp, 180._dp, 180._dp /)
INTEGER :: i
penalty = 1._dp
DO i = 1, npar
IF(pars_(i) < par_min(i)) THEN
pars_(i) = par_min(i)
penalty=penalty*10
ENDIF
IF(pars_(i) > par_max(i)) THEN
pars_(i) = par_max(i)
penalty=penalty*10
ENDIF
ENDDO
END SUBROUTINE
function rotx (theta) RESULT(R)
IMPLICIT NONE
REAL(DP),INTENT(in) :: theta
REAL(DP) :: R(3,3)
R(:,1) = (/ 1._dp, 0._dp, 0._dp/)
R(:,2) = (/ 0._dp, cos(theta), -sin(theta) /)
R(:,3) = (/ 0._dp, sin(theta), cos(theta) /)
endfunction
function roty (theta) RESULT(R)
IMPLICIT NONE
REAL(DP),INTENT(in) :: theta
REAL(DP) :: R(3,3)
R(:,1) = (/ cos(theta), 0._dp, sin(theta) /)
R(:,2) = (/ 0._dp,1._dp,0._dp/)
R(:,3) = (/ -sin(theta), 0._dp, cos(theta)/)
endfunction
function rotz (theta) RESULT(R)
IMPLICIT NONE
REAL(DP),INTENT(in) :: theta
REAL(DP) :: R(3,3)
R(:,1) = (/cos(theta),-sin(theta), 0._dp /)
R(:,2) = (/sin(theta), cos(theta), 0._dp /)
R(:,3) = (/0._dp,0._dp,1._dp /)
endfunction
function rot(alpha,beta,gamma) RESULT(R)
IMPLICIT NONE
REAL(DP),INTENT(in) :: alpha,beta,gamma
REAL(DP) :: R(3,3)
R = matmul(matmul(rotx(alpha),roty(beta)), rotz(gamma))
endfunction
END PROGRAM ibrav2cell
!----------------------------------------------------------------------
|
