← Back to branch summary

~siesta-ts/siesta/trunk_ts_soc

« back to all changes in this revision

Viewing changes to Src/diag3k.F

  • Committer: Nils Wittemeier
  • Date: 2019-02-14 07:45:07 UTC
  • mfrom: (746.1.15 trunk)
  • Revision ID: nils@4wittemeier.de-20190214074507-1mvzbmj9kw19gllr
Merged trunk 761

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/diag3k.F
179
179
 
180
180
      do ik = 1,nk
181
181
      
182
 
        Saux = dcmplx(0.0_dp,0.0_dp)
183
 
        Haux = dcmplx(0.0_dp,0.0_dp)
 
182
        Saux = cmplx(0.0_dp,0.0_dp,dp)
 
183
        Haux = cmplx(0.0_dp,0.0_dp,dp)
184
184
      
185
185
        do iuo = 1,nuo
186
186
          do j = 1,numh(iuo)
190
190
            kxij = kpoint(1,ik) * xij(1,ind) +
191
191
     .             kpoint(2,ik) * xij(2,ind) +
192
192
     .             kpoint(3,ik) * xij(3,ind)
193
 
            kphs = cdexp(dcmplx(0.0_dp, -1.0_dp)*kxij)
 
193
            kphs = exp(cmplx(0.0_dp, -1.0_dp, dp) * kxij)
194
194
      
195
195
            Saux(1,juo,1,iuo) = Saux(1,juo,1,iuo) + S(ind)   * kphs
196
196
            Saux(2,juo,2,iuo) = Saux(2,juo,2,iuo) + S(ind)   * kphs
197
197
            Haux(1,juo,1,iuo) = Haux(1,juo,1,iuo) + 
198
 
     .                            dcmplx(H(ind,1), H(ind,5)) * kphs
 
198
     .                            cmplx(H(ind,1), H(ind,5),dp) * kphs
199
199
            Haux(2,juo,2,iuo) = Haux(2,juo,2,iuo) + 
200
 
     .                            dcmplx(H(ind,2), H(ind,6)) * kphs
 
200
     .                            cmplx(H(ind,2), H(ind,6),dp) * kphs
201
201
            Haux(1,juo,2,iuo) = Haux(1,juo,2,iuo)
202
 
     .                        + dcmplx(H(ind,3), - H(ind,4)) * kphs
 
202
     .                        + cmplx(H(ind,3), - H(ind,4),dp) * kphs
203
203
            Haux(2,juo,1,iuo) = Haux(2,juo,1,iuo)
204
 
     .                        + dcmplx(H(ind,7), + H(ind,8)) * kphs
 
204
     .                        + cmplx(H(ind,7), + H(ind,8),dp) * kphs
205
205
          enddo
206
206
        enddo
207
207
      
252
252
        enddo
253
253
      
254
254
!       Find eigenvectors 
255
 
        Saux = dcmplx(0.0_dp,0.0_dp)
256
 
        Haux = dcmplx(0.0_dp,0.0_dp)
 
255
        Saux = cmplx(0.0_dp,0.0_dp,dp)
 
256
        Haux = cmplx(0.0_dp,0.0_dp,dp)
257
257
        do iuo = 1,nuo
258
258
          do j = 1,numh(iuo)
259
259
            ind = listhptr(iuo) + j
262
262
            kxij = kpoint(1,ik) * xij(1,ind) +
263
263
     .             kpoint(2,ik) * xij(2,ind) +
264
264
     .             kpoint(3,ik) * xij(3,ind)
265
 
            kphs = cdexp(dcmplx(0.0_dp, -1.0_dp)*kxij)
 
265
            kphs = exp(cmplx(0.0_dp, -1.0_dp, dp) * kxij)
266
266
      
267
267
            Saux(1,juo,1,iuo) = Saux(1,juo,1,iuo) + S(ind)   * kphs
268
268
            Saux(2,juo,2,iuo) = Saux(2,juo,2,iuo) + S(ind)   * kphs
269
269
            Haux(1,juo,1,iuo) = Haux(1,juo,1,iuo) + 
270
 
     .                            dcmplx(H(ind,1), H(ind,5)) * kphs
 
270
     .                            cmplx(H(ind,1), H(ind,5),dp) * kphs
271
271
            Haux(2,juo,2,iuo) = Haux(2,juo,2,iuo) + 
272
 
     .                            dcmplx(H(ind,2), H(ind,6)) * kphs
 
272
     .                            cmplx(H(ind,2), H(ind,6),dp) * kphs
273
273
            Haux(1,juo,2,iuo) = Haux(1,juo,2,iuo)
274
 
     .                        + dcmplx(H(ind,3), - H(ind,4)) * kphs
 
274
     .                        + cmplx(H(ind,3), - H(ind,4),dp) * kphs
275
275
            Haux(2,juo,1,iuo) = Haux(2,juo,1,iuo)
276
 
     .                        + dcmplx(H(ind,7), + H(ind,8)) * kphs
 
276
     .                        + cmplx(H(ind,7), + H(ind,8),dp) * kphs
277
277
          enddo
278
278
        enddo
279
279
 
286
286
          call die('Terminating due to failed diagonalisation')
287
287
        elseif (ierror.lt.0) then
288
288
!         Repeat diagonalisation with increased memory to handle clustering
289
 
          Saux = dcmplx(0.0_dp,0.0_dp)
290
 
          Haux = dcmplx(0.0_dp,0.0_dp)
 
289
          Saux = cmplx(0.0_dp,0.0_dp,dp)
 
290
          Haux = cmplx(0.0_dp,0.0_dp,dp)
291
291
          do iuo = 1,nuo
292
292
            do j = 1,numh(iuo)
293
293
              ind = listhptr(iuo) + j
296
296
              kxij = kpoint(1,ik) * xij(1,ind) +
297
297
     .               kpoint(2,ik) * xij(2,ind) +
298
298
     .               kpoint(3,ik) * xij(3,ind)
299
 
              kphs = cdexp(dcmplx(0.0_dp, -1.0_dp)*kxij)
 
299
              kphs = exp(cmplx(0.0_dp, -1.0_dp, dp) * kxij)
300
300
      
301
301
              Saux(1,juo,1,iuo) = Saux(1,juo,1,iuo) + S(ind)   * kphs
302
302
              Saux(2,juo,2,iuo) = Saux(2,juo,2,iuo) + S(ind)   * kphs
303
303
              Haux(1,juo,1,iuo) = Haux(1,juo,1,iuo) + 
304
 
     .                              dcmplx(H(ind,1), H(ind,5)) * kphs
 
304
     .                              cmplx(H(ind,1), H(ind,5),dp) * kphs
305
305
              Haux(2,juo,2,iuo) = Haux(2,juo,2,iuo) + 
306
 
     .                              dcmplx(H(ind,2), H(ind,6)) * kphs
 
306
     .                              cmplx(H(ind,2), H(ind,6),dp) * kphs
307
307
              Haux(1,juo,2,iuo) = Haux(1,juo,2,iuo)
308
 
     .                          + dcmplx(H(ind,3), - H(ind,4)) * kphs
 
308
     .                          + cmplx(H(ind,3), - H(ind,4),dp) * kphs
309
309
              Haux(2,juo,1,iuo) = Haux(2,juo,1,iuo)
310
 
     .                          + dcmplx(H(ind,7), + H(ind,8)) * kphs
 
310
     .                          + cmplx(H(ind,7), + H(ind,8),dp) * kphs
311
311
            enddo
312
312
          enddo
313
313
          call cdiag(Haux,Saux,2*nuotot,2*nuo,2*nuotot,caux,psi,
336
336
! D_{i,j}(1,2) = 0.5 ( D_{i,j}(1,2) + D_{i,j}(2,1)^* )
337
337
! can not be enforced here.
338
338
 
339
 
        Dk  = dcmplx(0.0_dp,0.0_dp)
340
 
        Ek  = dcmplx(0.0_dp,0.0_dp)
 
339
        Dk  = cmplx(0.0_dp,0.0_dp,dp)
 
340
        Ek  = cmplx(0.0_dp,0.0_dp,dp)
341
341
      
342
342
        BNode = 0
343
343
        iie = 0
347
347
            iie = iie + 1
348
348
          endif
349
349
      
350
 
          caux(:,:) = dcmplx(0.0_dp,0.0_dp)
 
350
          caux(:,:) = cmplx(0.0_dp,0.0_dp,dp)
351
351
          if (qe.gt.occtol) then
352
352
            if (Node.eq.BNode) then
353
353
              do j = 1,nuotot
365
365
              do juo = 1,nuotot
366
366
      
367
367
!------- 1,1 -----------------------------------------------------------
368
 
               cicj = caux(1,iio) * dconjg(caux(1,juo))
 
368
               cicj = caux(1,iio) * conjg(caux(1,juo))
369
369
               Dk(1,juo,1,iuo) = Dk(1,juo,1,iuo) + qe * cicj
370
370
               Ek(1,juo,1,iuo) = Ek(1,juo,1,iuo) + ee * cicj
371
371
!------- 2,2 -----------------------------------------------------------
372
 
               cicj = caux(2,iio) * dconjg(caux(2,juo))
 
372
               cicj = caux(2,iio) * conjg(caux(2,juo))
373
373
               Dk(2,juo,2,iuo) = Dk(2,juo,2,iuo) + qe * cicj
374
374
               Ek(2,juo,2,iuo) = Ek(2,juo,2,iuo) + ee * cicj
375
375
!------- 1,2 -----------------------------------------------------------
376
 
               cicj = caux(1,iio) * dconjg(caux(2,juo))
 
376
               cicj = caux(1,iio) * conjg(caux(2,juo))
377
377
               Dk(1,juo,2,iuo) = Dk(1,juo,2,iuo) + qe * cicj
378
378
               Ek(1,juo,2,iuo) = Ek(1,juo,2,iuo) + ee * cicj
379
379
!------- 2,1 -----------------------------------------------------------
380
 
               cicj = caux(2,iio) * dconjg(caux(1,juo))
 
380
               cicj = caux(2,iio) * conjg(caux(1,juo))
381
381
               Dk(2,juo,1,iuo) = Dk(2,juo,1,iuo) + qe * cicj
382
382
               Ek(2,juo,1,iuo) = Ek(2,juo,1,iuo) + ee * cicj
383
383
              enddo
398
398
          kxij = kpoint(1,ik) * xij(1,ind) +
399
399
     .           kpoint(2,ik) * xij(2,ind) +
400
400
     .           kpoint(3,ik) * xij(3,ind)
401
 
          kphs = cdexp(dcmplx(0.0_dp,-1.0_dp)*kxij)
 
401
          kphs = exp(cmplx(0.0_dp,-1.0_dp,dp) * kxij)
402
402
      
403
403
          D11 = Dk(1,juo,1,iuo) * kphs
404
404
          D22 = Dk(2,juo,2,iuo) * kphs
405
405
          D12 = Dk(1,juo,2,iuo) * kphs
406
406
          D21 = Dk(2,juo,1,iuo) * kphs
407
407
      
408
 
          Dnew(ind,1) = Dnew(ind,1) + dreal(D11)
409
 
          Dnew(ind,2) = Dnew(ind,2) + dreal(D22)
410
 
          Dnew(ind,3) = Dnew(ind,3) + dreal(D12)
411
 
          Dnew(ind,4) = Dnew(ind,4) - dimag(D12)
412
 
          Dnew(ind,5) = Dnew(ind,5) + dimag(D11) 
413
 
          Dnew(ind,6) = Dnew(ind,6) + dimag(D22)
414
 
          Dnew(ind,7) = Dnew(ind,7) + dreal(D21) 
415
 
          Dnew(ind,8) = Dnew(ind,8) + dimag(D21)
 
408
          Dnew(ind,1) = Dnew(ind,1) + real(D11,dp)
 
409
          Dnew(ind,2) = Dnew(ind,2) + real(D22,dp)
 
410
          Dnew(ind,3) = Dnew(ind,3) + real(D12,dp)
 
411
          Dnew(ind,4) = Dnew(ind,4) - aimag(D12)
 
412
          Dnew(ind,5) = Dnew(ind,5) + aimag(D11) 
 
413
          Dnew(ind,6) = Dnew(ind,6) + aimag(D22)
 
414
          Dnew(ind,7) = Dnew(ind,7) + real(D21,dp) 
 
415
          Dnew(ind,8) = Dnew(ind,8) + aimag(D21)
416
416
 
417
417
      
418
418
          D11 = Ek(1,juo,1,iuo) * kphs
420
420
          D12 = Ek(1,juo,2,iuo) * kphs
421
421
          D21 = Ek(2,juo,1,iuo) * kphs
422
422
      
423
 
          Enew(ind,1) = Enew(ind,1) + dreal(D11)
424
 
          Enew(ind,2) = Enew(ind,2) + dreal(D22)
425
 
          Enew(ind,3) = Enew(ind,3) + dreal(D12) 
426
 
          Enew(ind,4) = Enew(ind,4) - dimag(D12) 
 
423
          Enew(ind,1) = Enew(ind,1) + real(D11,dp)
 
424
          Enew(ind,2) = Enew(ind,2) + real(D22,dp)
 
425
          Enew(ind,3) = Enew(ind,3) + real(D12,dp) 
 
426
          Enew(ind,4) = Enew(ind,4) - aimag(D12) 
427
427
      
428
428
        enddo
429
429
        enddo