scnc
/
conquest
mirror of https://github.com/OrderN/CONQUEST-release.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

First pass for output levels 0-2 

Mostly tidied output for iprint levels up to 2, for
SCF, LFD and basis optimisation.  Also atomic
movement to some extent.
This commit is contained in:
David Bowler 2020-07-24 14:39:18 +01:00
parent 0c284534d4
commit b624857d23
23 changed files with 782 additions and 990 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/DMMinModule.f90
View File

@ -179,7 +179,8 @@ contains
ne_in_cell, ne_spin_in_cell, flag_dump_L, &
flag_SkipEarlyDM, flag_XLBOMD, &
flag_propagateX, flag_dissipation, &
integratorXL, runtype, flag_exx, flag_diagonalisation
integratorXL, runtype, flag_exx, &
flag_diagonalisation, min_layer
use mult_module, only: matrix_transpose, matT, matTtran, matL, &
matS, matrix_sum, matK
use McWeeny, only: InitMcW, McWMin

174
src/DiagModule.f90
View File

@ -143,7 +143,7 @@
module DiagModule
use datatypes
use global_module, only: io_lun, area_DM, iprint_DM, flag_diagonalisation
use global_module, only: io_lun, area_DM, iprint_DM, min_layer, flag_diagonalisation
use GenComms, only: cq_abort, inode, ionode, myid
use timer_module, only: start_timer, stop_timer
use timer_stdclocks_module, only: tmr_std_matrices
@ -475,7 +475,7 @@ contains
use datatypes
use numbers
use units
use global_module, only: iprint_DM, ni_in_cell, numprocs, &
use global_module, only: iprint_DM, min_layer, ni_in_cell, numprocs, &
area_DM, flag_fix_spin_population, &
nspin, spin_factor, flag_DeltaSCF, flag_excite, &
flag_local_excitation, dscf_HOMO_thresh, &
@ -532,7 +532,8 @@ contains
real(double), dimension(:),allocatable :: abs_wf
if (iprint_DM >= 2 .AND. myid == 0) &
min_layer = min_layer - 1
if (iprint_DM + min_layer >= 2 .AND. myid == 0) &
write (io_lun, fmt='(10x,"Entering FindEvals")')
prim_size = 0
do i = 1, bundle%n_prim
@ -550,7 +551,7 @@ contains
! First diagonalisation - get eigenvalues only (so that we can find Efermi)
time0 = mtime ()
if (iprint_DM >= 2 .and. (inode == ionode)) &
if (iprint_DM + min_layer >= 2 .and. (inode == ionode)) &
write (io_lun, fmt='(10x,"In FindEvals, tolerance is ", g20.12)') &
abstol
@ -571,7 +572,7 @@ contains
end do ! spin
! Allocate matrices to store band K matrices
time1 = mtime()
if (iprint_DM >= 2 .AND. myid == 0) &
if (iprint_DM + min_layer >= 2 .AND. myid == 0) &
write (io_lun, 2) myid, time1 - time0
! If we are trying to localise a level, we do NOT want to excite just yet
if(flag_DeltaSCF.AND.flag_excite.AND.flag_local_excitation) then
@ -619,7 +620,7 @@ contains
if(E_DOS_min>w(1,i,1)) E_DOS_min = w(1,i,1)
if(E_DOS_max<w(matrix_size,i,1)) E_DOS_max = w(matrix_size,i,1)
end do
if(myid==0.AND.iprint_DM>=2) write(io_lun,fmt='(2x,"DOS limits set automatically: ",2f12.5)') &
if(myid==0.AND.iprint_DM + min_layer>=2) write(io_lun,fmt='(2x,"DOS limits set automatically: ",2f12.5)') &
E_DOS_min, E_DOS_max
end if
dE_DOS = (E_DOS_max - E_DOS_min)/real(n_DOS-1,double)
@ -688,7 +689,7 @@ contains
end do
call gsum(setA)
call gsum(setB)
if(inode==ionode.AND.iprint_DM>2) &
if(inode==ionode.AND.iprint_DM + min_layer>2) &
write(io_lun,fmt='(4x,"DeltaSCF HOMO search: Band, coefficients A/B: ",i5,2f12.5)') band,setA,setB
if(setA>dscf_HOMO_thresh) then
dscf_source_level = band
@ -741,7 +742,7 @@ contains
end do
call gsum(setA)
call gsum(setB)
if(inode==ionode.AND.iprint_DM>2) &
if(inode==ionode.AND.iprint_DM + min_layer>2) &
write(io_lun,fmt='(4x,"DeltaSCF LUMO search: Band, coefficients A/B: ",i5,2f12.5)') band,setA,setB
if(setA>dscf_LUMO_thresh) then
dscf_target_level = band
@ -752,7 +753,7 @@ contains
call findFermi(electrons, w, matrix_size, nkp, Efermi, occ)
end if ! DeltaSCF localised excitation
! Now write out eigenvalues and occupancies
if (iprint_DM == 2 .AND. myid == 0) then
if (iprint_DM + min_layer == 2 .AND. myid == 0) then
bandE = zero
do i = 1, nkp
write (io_lun, 7) i, kk(1,i), kk(2,i), kk(3,i)
@ -790,7 +791,7 @@ contains
write(io_lun, 4) en_conv * two * bandE(1), en_units(energy_units)
end if
end do ! do i = 1, nkp
else if (iprint_DM >= 3 .AND. myid == 0) then
else if (iprint_DM + min_layer >= 3 .AND. myid == 0) then
bandE = zero
do i = 1, nkp
write (io_lun, 7) i, kk(1,i), kk(2,i), kk(3,i)
@ -814,7 +815,7 @@ contains
write(io_lun, 4) en_conv * two * bandE(1), en_units(energy_units)
end if
end do ! do i = 1, nkp
end if ! if(iprint_DM>=1.AND.myid==0)
end if ! if(iprint_DM + min_layer>=1.AND.myid==0)
time0 = mtime()
do spin = 1, nspin
@ -869,7 +870,7 @@ contains
end if
end if
! Build K and K_dn from the eigenvectors
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, *) myid, ' Calling buildK ', &
Hrange, matK(spin)
! Output wavefunction coefficients
@ -894,7 +895,7 @@ contains
entropy_local(spin) = &
locc(spin) * log(locc(spin)) + &
(one - locc(spin)) * log (one - locc(spin))
if (iprint_DM > 3 .and. inode == ionode) &
if (iprint_DM + min_layer > 3 .and. inode == ionode) &
write (io_lun, &
fmt='(2x,"Spin, Occ, wt: ", i1, 2f12.8, &
&" ent: ", f20.12)') &
@ -962,26 +963,26 @@ contains
deallocate(total_DOS)
end if
if (iprint_DM > 3 .and. inode == ionode) &
if (iprint_DM + min_layer > 3 .and. inode == ionode) &
write (io_lun, *) "Entropy, TS: ", entropy, kT * entropy
! store entropy as TS instead of S
entropy = entropy * kT
time1 = mtime()
if (iprint_DM >= 2 .and. inode == ionode) &
if (iprint_DM + min_layer >= 2 .and. inode == ionode) &
write (io_lun, 3) myid, time1 - time0
! -------------------------------------------------------------
! End diagonalisation
! -------------------------------------------------------------
! Write out the Fermi Energy
if (iprint_DM >= 1 .and. inode == ionode) then
if (iprint_DM + min_layer >= 1 .and. inode == ionode) then
do spin = 1, nspin
write (io_lun, 13) spin, en_conv * Efermi(spin), &
en_units(energy_units)
end do
end if
! Write out the band energy and trace of K
if (iprint_DM >= 1) then
if (iprint_DM + min_layer >= 1) then
! for tr(K.H)
bandE_total = zero
do spin = 1, nspin
@ -1025,7 +1026,7 @@ contains
call reg_dealloc_mem(area_DM, matrix_size * prim_size * nspin, type_cplx)
! global
call endDiag
min_layer = min_layer + 1
return
2 format(10x,'Proc: ',i5, ' Time taken for eval diag: ',f20.8,' ms')
@ -2387,7 +2388,7 @@ contains
use datatypes
use numbers
use global_module, only: iprint_DM, nspin
use global_module, only: iprint_DM, nspin, min_layer
use GenComms, only: myid
implicit none
@ -2413,7 +2414,7 @@ contains
electrons_total = two * electrons(1)
end if
if (iprint_DM >= 2 .and. (inode == ionode)) then
if (iprint_DM + min_layer >= 2 .and. (inode == ionode)) then
if (nspin == 1) then
write (io_lun, 1) myid, electrons_total
else
@ -2424,7 +2425,7 @@ contains
! find the correct bracket trapping Ef
labspin: do spin = 1, nspin
if (iprint_DM >= 2 .and. nspin == 2 .and. inode == ionode) &
if (iprint_DM + min_layer >= 2 .and. nspin == 2 .and. inode == ionode) &
write (io_lun, 3) myid, spin
select case (flag_smear_type)
@ -2444,7 +2445,7 @@ contains
if (thisElec(spin) < electrons(spin)) then ! found a lower bound
if (iprint_DM >= 4 .and. (inode == ionode)) &
if (iprint_DM + min_layer >= 4 .and. (inode == ionode)) &
write (io_lun, 4) myid, Ef(spin)
lowEf(spin) = Ef(spin)
lowElec(spin) = thisElec(spin)
@ -2465,14 +2466,14 @@ contains
lowElec(spin) = highElec(spin)
highEf(spin) = highEf(spin) + incEf(spin)
call occupy(occ, eig, highEf, highElec, nbands, nkp, spin=spin)
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 5) myid, highEf(spin), highElec(spin)
ibrkt = ibrkt + 1
end do ! while (highElec(spin) < electrons(spin))
else ! found an upper bound
if (iprint_DM >= 4 .AND. inode == ionode) &
if (iprint_DM + min_layer >= 4 .AND. inode == ionode) &
write (io_lun, 6) myid, Ef(spin)
highEf(spin) = Ef(spin)
highElec(spin) = thisElec(spin)
@ -2493,14 +2494,14 @@ contains
highElec(spin) = lowElec(spin)
lowEf(spin) = lowEf(spin) - incEf(spin)
call occupy(occ, eig, lowEf, lowElec, nbands, nkp, spin=spin)
if (iprint_DM >= 4 .AND. inode == ionode) &
if (iprint_DM + min_layer >= 4 .AND. inode == ionode) &
write (io_lun, 5) myid, lowEf(spin), lowElec(spin)
ibrkt = ibrkt + 1
end do
end if ! if (thisElec(spin) < electrons(spin))
if (iprint_DM > 3 .and. inode == ionode) &
if (iprint_DM + min_layer > 3 .and. inode == ionode) &
write (io_lun, 12) myid, lowEf(spin), highEf(spin)
case (1) ! Methfessel-Paxton smearing
@ -2563,11 +2564,11 @@ contains
lowElec(spin) = highElec(spin)
highEf(spin) = lowEf(spin) + incEf(spin)
call occupy(occ, eig, highEf, highElec, nbands, nkp, spin=spin)
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 5) myid, highEf(spin), highElec(spin)
end do
if (iprint_DM > 3 .and. inode == ionode) &
if (iprint_DM + min_layer > 3 .and. inode == ionode) &
write (io_lun, 12) myid, lowEf(spin), highEf(spin)
case default
@ -2593,7 +2594,7 @@ contains
call occupy(occ, eig, Ef, thisElec, nbands, nkp, spin=spin)
end do
if (iprint_DM >= 2 .AND. inode == ionode) then
if (iprint_DM + min_layer >= 2 .AND. inode == ionode) then
if (nspin == 1) then
write (io_lun, 10) Ef(spin)
else
@ -2659,7 +2660,7 @@ contains
use datatypes
use numbers
use global_module, only: iprint_DM, nspin, spin_factor
use global_module, only: iprint_DM, nspin, spin_factor, min_layer
use GenComms, only: myid
implicit none
@ -2683,7 +2684,7 @@ contains
case (0) ! Fermi smearing
if (iprint_DM >= 2 .AND. inode == ionode) &
if (iprint_DM + min_layer >= 2 .AND. inode == ionode) &
write (io_lun, 1) myid, electrons_total
! Take first guess as double filling each band at first k point
ne = int(electrons_total / two)
@ -2698,7 +2699,7 @@ contains
! Find two values than bracket true Ef
incEf(:) = one
if (thisElec < electrons_total) then ! found lower bound
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 2) myid, Ef(1) ! Ef(1) = Ef(2) always
lowEf(:) = Ef(:)
lowElec = thisElec
@ -2719,12 +2720,12 @@ contains
highEf(:) = highEf(:) + incEf(:)
call occupy(occ, eig, highEf, electrons, nbands, nkp)
highElec = spin_factor * sum(electrons(:))
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 3) myid, highEf(1), highElec
ibrkt = ibrkt + 1
end do
else ! found upper bound
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 6) myid, Ef(1)
highEf(:) = Ef(:)
highElec = thisElec
@ -2745,12 +2746,12 @@ contains
lowEf(:) = lowEf(:) - incEf(:)
call occupy(occ, eig, lowEf, electrons, nbands, nkp)
lowElec = spin_factor * sum(electrons(:))
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 3) myid, lowEf(1), lowElec
ibrkt = ibrkt + 1
end do
end if
if (iprint_DM > 3 .AND. inode == ionode) &
if (iprint_DM + min_layer > 3 .AND. inode == ionode) &
write (io_lun, 5) myid, lowEf(1), highEf(1)
! search method to use in the case of Methmessel Paxton smearing
@ -2816,7 +2817,7 @@ contains
highEf(:) = lowEf(:) + incEf(:)
call occupy(occ, eig, highEf, electrons, nbands, nkp)
highElec = spin_factor * sum(electrons(:))
if (iprint_DM >= 4 .and. inode == ionode) &
if (iprint_DM + min_layer >= 4 .and. inode == ionode) &
write (io_lun, 3) myid, highEf(1), highElec
end do
@ -2845,7 +2846,7 @@ contains
thisElec = spin_factor * sum(electrons(:))
end do
if (iprint_DM >= 2 .and. inode == ionode) write (io_lun, 8) Ef(1)
if (iprint_DM + min_layer >= 2 .and. inode == ionode) write (io_lun, 8) Ef(1)
return
@ -2916,7 +2917,7 @@ contains
subroutine occupy(occu, ebands, Ef, electrons, nbands, nkp, spin)
use datatypes
use numbers
use global_module, only: iprint_DM, nspin
use global_module, only: iprint_DM, nspin, min_layer
use GenComms, only: myid
implicit none
@ -2962,7 +2963,7 @@ contains
end do kp
end do labspin
if (iprint_DM >=5 .and. (inode == ionode)) then
if (iprint_DM + min_layer >=5 .and. (inode == ionode)) then
if (nspin == 1) then
write (io_lun, 1) myid, Ef(1), two * electrons(1)
else if (present (spin)) then
@ -3290,7 +3291,7 @@ contains
matrix_size
use global_module, only: numprocs, iprint_DM, id_glob, &
ni_in_cell, x_atom_cell, y_atom_cell, &
z_atom_cell, max_wf
z_atom_cell, max_wf, min_layer
use mpi
use GenBlas, only: dot
use GenComms, only: myid
@ -3335,7 +3336,7 @@ contains
real(double) :: occ_correction
call start_timer(tmr_std_matrices)
if(iprint_DM>3.AND.myid==0) write(io_lun,fmt='(10x,"Entering &
if(iprint_DM + min_layer>3.AND.myid==0) write(io_lun,fmt='(10x,"Entering &
&buildK ",i4)') matA
! get occ_correction
@ -3355,17 +3356,17 @@ contains
send_prim = 0
num_send = 0
norb_send = 0
if(iprint_DM>3.AND.myid==0) write(io_lun,*) 'buildK: Stage one'
if(iprint_DM + min_layer>3.AND.myid==0) write(io_lun,*) 'buildK: Stage one'
! Step one - work out which processors we need to exchange data with
do part = 1,bundle%groups_on_node ! Loop over primary set partitions
if(iprint_DM>=5.AND.myid==0) write(io_lun,1) myid,part
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,1) myid,part
if(bundle%nm_nodgroup(part)>0) then ! If there are atoms in partition
CC = parts%ngnode(parts%inode_beg(myid+1)+part-1)
do memb = 1,bundle%nm_nodgroup(part) ! Loop over atoms
if(iprint_DM>=5.AND.myid==0) write(io_lun,2) myid,memb
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,2) myid,memb
prim_atom = bundle%nm_nodbeg(part)+memb-1
do neigh = 1, mat(part,range)%n_nab(memb) ! Loop over neighbours of atom
if(iprint_DM>=5.AND.myid==0) write(io_lun,3) myid,neigh
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,3) myid,neigh
ist = mat(part,range)%i_acc(memb)+neigh-1
! Establish FSC number of neighbour
Col_FSC_part = BCS_parts%lab_cell(mat(part,range)%i_part(ist))
@ -3376,12 +3377,12 @@ contains
FSC_atom = id_glob(parts%icell_beg(Col_FSC_part)+Col_FSC_seq-1)
! Find if we have seen this before
flag = .false.
if(iprint_DM>=5.AND.myid==0) write(io_lun,*) 'prim, neigh, FSC: ',prim_atom, neigh, FSC_atom
if(iprint_DM>=5.AND.myid==0) write(io_lun,*) 'curr_loc_atoms: ',current_loc_atoms(owning_proc)
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,*) 'prim, neigh, FSC: ',prim_atom, neigh, FSC_atom
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,*) 'curr_loc_atoms: ',current_loc_atoms(owning_proc)
if(current_loc_atoms(owning_proc)>0) then
do i=1,current_loc_atoms(owning_proc)
if(atom_list(owning_proc,i)==FSC_atom) then
if(iprint_DM>=5.AND.myid==0) write(io_lun,*) 'Loc atom: ',i, LocalAtom(FSC_atom)
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,*) 'Loc atom: ',i, LocalAtom(FSC_atom)
ints(owning_proc,LocalAtom(FSC_atom)) = ints(owning_proc,LocalAtom(FSC_atom)) + 1
send_prim(owning_proc,prim_atom) = nsf_species(bundle%species(prim_atom))
flag = .true.
@ -3403,22 +3404,22 @@ contains
end if ! End if nm_nodgroup > 0
end do ! End do part=1,groups_on_node
! Find max value of current_loc_atoms and interactions
if(iprint_DM>3.AND.myid==0) write(io_lun,*) 'buildK: Stage two'
if(iprint_DM + min_layer>3.AND.myid==0) write(io_lun,*) 'buildK: Stage two'
maxloc = 0
maxint = 0
maxsend = 0
do i=1,numprocs
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) myid,' Curr loc atoms: ',i,current_loc_atoms(i)
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) myid,' Curr loc atoms: ',i,current_loc_atoms(i)
if(current_loc_atoms(i)>maxloc) maxloc = current_loc_atoms(i)
do j=1,bundle%mx_iprim ! Needs to be mx_iprim because goes over primary atoms on REMOTE processors
if(ints(i,j)>maxint) maxint = ints(i,j)
if(send_prim(i,j)>0) num_send(i) = num_send(i) + 1
norb_send(i) = norb_send(i) + send_prim(i,j)
if(iprint_DM>=5.AND.myid==0) write(io_lun,4) myid,j,send_prim(i,j),num_send(i)
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,4) myid,j,send_prim(i,j),num_send(i)
end do
if(num_send(i)>maxsend) maxsend = num_send(i)
end do
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) myid,' Maxima: ',maxloc, maxint, maxsend
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) myid,' Maxima: ',maxloc, maxint, maxsend
! Allocate recv_info
allocate(send_info(numprocs,maxsend),send_orbs(numprocs,maxsend),send_off(numprocs,maxsend), &
prim_orbs(bundle%mx_iprim),STAT=stat)
@ -3462,14 +3463,14 @@ contains
if(stat/=0) call cq_abort('buildK: Error allocating recv_info !',stat)
end do
do part = 1,bundle%groups_on_node ! Loop over primary set partitions
if(iprint_DM>=5.AND.myid==0) write(io_lun,1) myid,part
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,1) myid,part
if(bundle%nm_nodgroup(part)>0) then ! If there are atoms in partition
CC = parts%ngnode(parts%inode_beg(myid+1)+part-1)
do memb = 1,bundle%nm_nodgroup(part) ! Loop over atoms
if(iprint_DM>=5.AND.myid==0) write(io_lun,2) myid,memb
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,2) myid,memb
prim_atom = bundle%nm_nodbeg(part)+memb-1
do neigh = 1, mat(part,range)%n_nab(memb) ! Loop over neighbours of atom
if(iprint_DM>=5.AND.myid==0) write(io_lun,3) myid,neigh
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,3) myid,neigh
ist = mat(part,range)%i_acc(memb)+neigh-1
! Establish FSC number of neighbour
Col_FSC_part = BCS_parts%lab_cell(mat(part,range)%i_part(ist))
@ -3480,10 +3481,10 @@ contains
FSC_atom = id_glob(parts%icell_beg(Col_FSC_part)+Col_FSC_seq-1)
! Work out a map from primary atom + FSC + identifier to distance and position in data_Matrix
locatom = LocalAtom(FSC_atom) ! Which atom in the list on the remote proc is this ?
if(iprint_DM>=5.AND.myid==0) write(io_lun,*) myid,' own, FSC, loc: ',owning_proc, FSC_atom, locatom, &
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,*) myid,' own, FSC, loc: ',owning_proc, FSC_atom, locatom, &
recv_info(owning_proc)%ints(locatom)
recv_info(owning_proc)%ints(locatom) = recv_info(owning_proc)%ints(locatom) + 1
if(iprint_DM>=5.AND.myid==0) write(io_lun,*) myid,' ints: ',recv_info(owning_proc)%ints(locatom)
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,*) myid,' ints: ',recv_info(owning_proc)%ints(locatom)
gcspart = BCS_parts%icover_ibeg(mat(part,range)%i_part(ist))+mat(part,range)%i_seq(ist)-1
!recv_info(owning_proc)%ndimi(locatom) = mat(part,range)%ndimi(memb)
recv_info(owning_proc)%ndimj(locatom) = mat(part,range)%ndimj(ist)
@ -3509,11 +3510,11 @@ contains
do i=1,matrix_size ! Effectively all bands
if(abs(occs(i))>RD_ERR) then
len = i !len+1
if(myid==0.AND.iprint_DM>=4) write(io_lun,*) 'Occ is ',occs(i)
if(myid==0.AND.iprint_DM + min_layer>=4) write(io_lun,*) 'Occ is ',occs(i)
end if
end do
len_occ = len
if(iprint_DM>3.AND.myid==0) write(io_lun,*) 'buildK: Stage three len:',len, matA
if(iprint_DM + min_layer>3.AND.myid==0) write(io_lun,*) 'buildK: Stage three len:',len, matA
! Step three - loop over processors, send and recv data and build K
allocate(send_fsc(bundle%mx_iprim),recv_to_FSC(bundle%mx_iprim),mapchunk(bundle%mx_iprim),STAT=stat)
if(stat/=0) call cq_abort('buildK: Error allocating send_fsc, recv_to_FSC and mapchunk',stat)
@ -3528,11 +3529,11 @@ contains
do i=1,numprocs
send_size = len*norb_send(send_proc+1)!num_send(send_proc+1)*nsf
recv_size = len*recv_info(recv_proc+1)%orbs!current_loc_atoms(recv_proc+1)*nsf
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Send and recv sizes: ',send_size, recv_size
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Send and recv sizes: ',send_size, recv_size
! Fill SendBuffer
allocate(SendBuffer(len,norb_send(send_proc+1)),STAT=stat)
if(stat/=0) call cq_abort('buildK: Unable to allocate SendBuffer !',stat)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Filling SendBuffer'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Filling SendBuffer'
orb_count = 0
do j=1,num_send(send_proc+1)
do nsf1=1,send_orbs(send_proc+1,j)
@ -3541,11 +3542,11 @@ contains
end do
! We also need to send a list of what FSC each primary atom sent corresponds to - use bundle%ig_prim
send_FSC(j) = bundle%ig_prim(send_info(send_proc+1,j))
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Building send_FSC: ',send_info(send_proc+1,j), &
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Building send_FSC: ',send_info(send_proc+1,j), &
bundle%ig_prim(send_info(send_proc+1,j)),send_FSC(j)
end do
if(orb_count/=norb_send(send_proc+1)) call cq_abort("Orbital mismatch in buildK: ",orb_count,norb_send(send_proc+1))
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Sending'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Sending'
! Now send
if(send_size>0) then
if(send_proc/=myid) then
@ -3554,36 +3555,37 @@ contains
end if
end if
! Now receive data
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Alloc RecvBuffer ',len,recv_info(recv_proc+1)%orbs
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Alloc RecvBuffer ',len,recv_info(recv_proc+1)%orbs
!allocate(RecvBuffer(len,current_loc_atoms(recv_proc+1)*nsf),STAT=stat)
allocate(RecvBuffer(len,recv_info(recv_proc+1)%orbs),STAT=stat)
if(stat/=0) call cq_abort('buildK: Unable to allocate RecvBuffer !',stat)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Recving'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Recving'
if(recv_size>0) then
if(recv_proc/=myid) then
call MPI_recv(recv_to_FSC,current_loc_atoms(recv_proc+1),MPI_INTEGER,recv_proc,recvtag,MPI_COMM_WORLD,mpi_stat,ierr)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Got recv_to_FSC'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Got recv_to_FSC'
call MPI_recv(RecvBuffer,recv_size,MPI_DOUBLE_COMPLEX,&
recv_proc,recvtag+1,MPI_COMM_WORLD,mpi_stat,ierr)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Got RecvBuffer'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Got RecvBuffer'
else
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'On-proc: getting recv_to_FSC'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'On-proc: getting recv_to_FSC'
recv_to_FSC(1:current_loc_atoms(recv_proc+1)) = send_FSC(1:current_loc_atoms(recv_proc+1))
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'On-proc: getting RecvBuffer'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'On-proc: getting RecvBuffer'
RecvBuffer(1:len,1:recv_info(recv_proc+1)%orbs) = SendBuffer(1:len,1:recv_info(recv_proc+1)%orbs)
end if
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Doing the mapchunk', recv_to_FSC
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Doing the mapchunk', recv_to_FSC
do j=1,current_loc_atoms(recv_proc+1)
mapchunk(j) = LocalAtom(recv_to_FSC(j))
end do
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'filling buffer'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'filling buffer'
do j=1,len_occ ! This is a loop over eigenstates
RecvBuffer(j,1:recv_info(recv_proc+1)%orbs) = RecvBuffer(j,1:recv_info(recv_proc+1)%orbs)*occ_correction*occs(j)
end do
orb_count = 0
do atom = 1,current_loc_atoms(recv_proc+1)
locatom = mapchunk(atom)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Atom, loc: ',atom,locatom,recv_info(recv_proc+1)%ints(locatom)
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Atom, loc: ',atom, &
locatom,recv_info(recv_proc+1)%ints(locatom)
! Scale the eigenvector coefficients we've received
! The factor of 0.5 is because the occupation numbers are from 0->2 (we expect 0->1 in K)
! The occupation numbers contain the k-point weight
@ -3596,10 +3598,12 @@ contains
! N.B. the routine used for dot is zdotc which takes the complex conjugate of the first vector
do inter = 1,recv_info(recv_proc+1)%ints(locatom)
prim = recv_info(recv_proc+1)%prim_atom(inter,locatom)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Inter: ',inter,prim
phase = kps(1)*recv_info(recv_proc+1)%dx(inter,locatom) + kps(2)*recv_info(recv_proc+1)%dy(inter,locatom) + &
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Inter: ',inter,prim
phase = kps(1)*recv_info(recv_proc+1)%dx(inter,locatom) + &
kps(2)*recv_info(recv_proc+1)%dy(inter,locatom) + &
kps(3)*recv_info(recv_proc+1)%dz(inter,locatom)
if(iprint_DM>=5.AND.myid==0) write(io_lun,*) 'Prim, where, phase: ',prim, whereMat, phase
if(iprint_DM + min_layer>=5.AND.myid==0) write(io_lun,*) 'Prim, where, phase: ', &
prim, whereMat, phase
rfac = cos(phase)
ifac = sin(phase)
do row_sup = 1,recv_info(recv_proc+1)%ndimj(locatom)
@ -3615,20 +3619,20 @@ contains
orb_count = orb_count + recv_info(recv_proc+1)%ndimj(locatom)
end do ! atom=current_loc_atoms
end if ! recv_size>0
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Calling MPI_Wait'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Calling MPI_Wait'
if(send_size>0.AND.myid/=send_proc) then
call MPI_Wait(req1,mpi_stat,ierr)
call MPI_Wait(req2,mpi_stat,ierr)
end if
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Calling dealloc'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Calling dealloc'
deallocate(RecvBuffer,STAT=stat)
if(stat/=0) call cq_abort("buildK: Failed to dealloc buffer",stat)
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Calling dealloc'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Calling dealloc'
deallocate(SendBuffer,STAT=stat)
if(stat/=0) call cq_abort("buildK: Failed to dealloc buffer",stat)
! Increment/decrement recv and send, and wrap
! Remember that we go from 0->numprocs-1
if(iprint_DM>=4.AND.myid==0) write(io_lun,*) 'Doing proc thang'
if(iprint_DM + min_layer>=4.AND.myid==0) write(io_lun,*) 'Doing proc thang'
send_proc = send_proc +1
if(send_proc.GT.numprocs-1) send_proc = 0
recv_proc = recv_proc -1
@ -4121,7 +4125,7 @@ contains
use datatypes
use numbers
use global_module, only: iprint_DM, flag_SpinDependentSF
use global_module, only: iprint_DM, flag_SpinDependentSF, min_layer
use mult_module, only: matH, matS
use ScalapackFormat, only: matrix_size, proc_rows, proc_cols, &
nkpoints_max, pgid, N_kpoints_in_pg, pg_kpoints, N_procs_in_pg, proc_groups
@ -4148,7 +4152,7 @@ contains
orfac = -one
il = 0
iu = 0
if (iprint_DM > 3 .and. (inode == ionode)) &
if (iprint_DM + min_layer > 3 .and. (inode == ionode)) &
write (io_lun, *) myid, ' Calling DistributeCQ_to_SC for H'
! Form the Hamiltonian and overlap for this k-point and send them to appropriate processors
if(PRESENT(kpassed)) then
@ -4160,7 +4164,7 @@ contains
call DistributeCQ_to_SC(DistribS, matS(spin_SF), index_kpoint, SCSmat(:,:,spin))
end if
! Now, if this processor is involved, do the diagonalisation
if (iprint_DM > 3 .and. inode == ionode) &
if (iprint_DM + min_layer > 3 .and. inode == ionode) &
write (io_lun, *) myid, 'Proc row, cols, me: ', &
proc_rows, proc_cols, me, index_kpoint, nkpoints_max
if (index_kpoint <= N_kpoints_in_pg(pgid)) then
@ -4263,14 +4267,14 @@ contains
! check if we get the correct map
call blacs_gridinfo(context, numrows, numcols, merow, mecol)
if (iprint_DM > 3 .AND. myid == 0) then
if (iprint_DM + min_layer > 3 .AND. myid == 0) then
write (io_lun, fmt="(10x, 'process_grid info: ', i5, i5)") numrows, numcols
write (io_lun, 1) myid, me, merow, mecol
end if
! Sizes of local "chunk", used to initialise submatrix info for ScaLAPACK
row_size = proc_start(myid+1)%rows * block_size_r
col_size = proc_start(myid+1)%cols * block_size_c
if (iprint_DM > 3 .AND. myid == 0) write (io_lun, 12) myid, row_size, col_size
if (iprint_DM + min_layer > 3 .AND. myid == 0) write (io_lun, 12) myid, row_size, col_size
! Register the description of the distribution of H
call descinit(desca, matrix_size, matrix_size, block_size_r, block_size_c, 0, 0, context, row_size, info)

40
src/H_matrix_module.f90
View File

@ -215,7 +215,7 @@ contains
area_ops, nspin, nspin_SF, &
spin_factor, blips, &
flag_analytic_blip_int, &
flag_neutral_atom
flag_neutral_atom, min_layer
use functions_on_grid, only: atomfns, H_on_atomfns, &
gridfunctions
use io_module, only: dump_matrix, dump_blips, &
@ -319,10 +319,10 @@ contains
!
! from here on, workspace support becomes h_on_atomfns...
! in fact, what we are getting here is (H_local - T) acting on support
call get_h_on_atomfns(iprint_ops, fixed_potential, electrons, rho, size)
call get_h_on_atomfns(iprint_ops + min_layer, fixed_potential, electrons, rho, size)
!
!
if (inode == ionode .and. iprint_ops > 2) &
if (inode == ionode .and. iprint_ops > 3) &
write (io_lun, *) 'Doing integration'
! Do the integration - support holds <phi| and workspace_support
! holds H|phi>. Inode starts from 1, and myid starts from 0.
@ -332,7 +332,7 @@ contains
matHatomf(spin), atomfns, &
H_on_atomfns(spin))
end do
if (inode == ionode .and. iprint_ops > 2) write (io_lun, *) 'Done integration'
if (inode == ionode .and. iprint_ops > 3) write (io_lun, *) 'Done integration'
!
!
if (iprint_ops > 4) then
@ -361,29 +361,29 @@ contains
!if (inode==ionode) print*, 'exx_pulay_r0 = ', exx_pulay_r0
if ( exx_niter < exx_siter ) then
! For first H building use pure DFT. To be improved for Hartree-Fock
if (inode == ionode .and. iprint_exx > 2) &
if (inode == ionode .and. iprint_exx > 3) &
write (io_lun, *) 'EXX: first guess from DFT'
!
else
!
if (inode == ionode .and. iprint_exx > 2) &
if (inode == ionode .and. iprint_exx > 3) &
write (io_lun, *) 'EXX: setting get_X_matrix'
call get_X_params(backtrace_level)
call exx_global_write()
!
!if (inode == ionode .and. iprint_exx > 2) &
!if (inode == ionode .and. iprint_exx > 3) &
!write (io_lun, *) 'EXX: doing get_X_matrix'
do spin = 1, nspin
if (inode == ionode .and. iprint_exx > 2) &
if (inode == ionode .and. iprint_exx > 3) &
write (io_lun, *) 'EXX: doing get_X_matrix: ', print_exxspin(spin)
call get_X_matrix(spin,backtrace_level)
end do
!
!if (inode == ionode .and. iprint_exx > 2) &
!if (inode == ionode .and. iprint_exx > 3) &
!write (io_lun, *) 'EXX: done get_X_matrix'
do spin = 1, nspin
if (inode == ionode .and. iprint_exx > 2) &
if (inode == ionode .and. iprint_exx > 3) &
write (io_lun, *) 'EXX: done get_X_matrix: ', print_exxspin(spin)
call matrix_sum(one, matHatomf(spin),-exx_alpha*half, matXatomf(spin))
end do
@ -443,6 +443,24 @@ contains
endif
!
!
! dump charges if required
if (flag_DumpChargeDensity .or. iprint_SC > 3) then
if(nspin==1) then
allocate(rho_total(size), STAT=stat)
if (stat /= 0) call cq_abort("Error allocating rho_total: ", size)
call reg_alloc_mem(area_ops, size, type_dbl)
rho_total(:) = spin_factor * rho(:,1)
call dump_charge(rho_total, size, inode, spin=0)
deallocate(rho_total, STAT=stat)
if (stat /= 0) call cq_abort("Error deallocating rho_total")
call reg_dealloc_mem(area_ops, size, type_dbl)
else if (nspin == 2) then
call dump_charge(rho(:,1), size, inode, spin=1)
call dump_charge(rho(:,2), size, inode, spin=2)
end if
end if
!
!
! Store the new H matrix for cDFT
if (flag_perform_cDFT) then
do spin = 1, nspin
@ -674,7 +692,7 @@ contains
end if
electrons_tot = spin_factor * sum(electrons(:))
if (inode == ionode .and. output_level >= 2) then
if (inode == ionode .and. output_level >= 3) then
write (io_lun, '(10x,a)') &
'get_h_on_atomfns: Electron Count, up, down and total:'
write (io_lun, '(10x, 3f25.15)') &

12
src/S_matrix_module.f90
View File

@ -134,7 +134,7 @@ contains
blips, PAOs, atomf, sf, &
ni_in_cell, nspin_SF, &
flag_perform_cdft, &
IPRINT_TIME_THRES1
IPRINT_TIME_THRES1, min_layer
use matrix_data, only: Srange
use mult_module, only: matS, matSatomf, AtomF_to_SF_transform
use io_module, only: dump_matrix
@ -199,7 +199,7 @@ contains
end if
! get the new InvS matrix
call Iter_Hott_InvS(iprint_ops, InvSMaxSteps, &
call Iter_Hott_InvS(iprint_ops + min_layer, InvSMaxSteps, &
InvSDeltaOmegaTolerance, ni_in_cell, &
inode, ionode, flag_do_SFtransform)
@ -747,7 +747,7 @@ contains
type(matrix_store_global) :: InfoGlob
if (atomf.ne.sf .and. .not.flag_do_SFtransform) then
if (inode.eq.ionode.and.output_level>=2) write(io_lun,*) &
if (inode.eq.ionode.and.output_level>=3) write(io_lun,*) &
'Now we have only Spao but not Ssf yet, so InvS is not calculated at present.'
else
matI = allocate_temp_matrix(TSrange,0)
@ -777,10 +777,10 @@ contains
n_orbs = n_orbs + real(nsf_species(species(i)),double)
end do
! First construct the identity
if (inode.eq.ionode.and.output_level>=2) write(io_lun,*) 'Zeroing data'
if (inode.eq.ionode.and.output_level>=3) write(io_lun,*) 'Zeroing data'
call matrix_scale(zero,matI)
call matrix_scale(zero,matT(1))
if (inode.eq.ionode.and.output_level>=2) write(io_lun,*) 'Creating I'
if (inode.eq.ionode.and.output_level>=3) write(io_lun,*) 'Creating I'
ip = 1
nb = 1
do np = 1,bundle%groups_on_node
@ -822,7 +822,7 @@ contains
enddo
call gsum(tot)
eps = 1.0_double/(tot)
if(output_level>1.and.inode==ionode) write(io_lun,*) 'Eps, tot: ',eps,tot
if(output_level>2.and.inode==ionode) write(io_lun,*) 'Eps, tot: ',eps,tot
call matrix_scale(zero,matT(spin_SF))
call matrix_sum(zero,matT(spin_SF),eps,matS(spin_SF))
enddo ! spin_SF

14
src/ScalapackFormat.f90
View File

@ -275,7 +275,7 @@ contains
if(myid==0.AND.iprint_DM>3) write(io_lun,1) blocks_r,blocks_c
maxrow = floor(real(blocks_r/proc_rows))+1
maxcol = floor(real(blocks_c/proc_cols))+1
if(iprint_DM>1.AND.myid==0) write(io_lun,*) 'maxrow, maxcol: ',maxrow,maxcol
if(iprint_DM>3.AND.myid==0) write(io_lun,*) 'maxrow, maxcol: ',maxrow,maxcol
call start_timer(tmr_std_allocation)
allocate(mapx(numprocs,maxrow,maxcol),mapy(numprocs,maxrow,maxcol),STAT=stat)
if(stat/=0) call cq_abort("ScalapackFormat: Could not alloc map",stat)
@ -492,9 +492,9 @@ contains
integer :: i, j, n, nrow, ncol, prow, pcol, proc
integer :: ng
if(iprint_DM>2.AND.myid==0) write(io_lun,*) myid,' Starting Ref To SC Blocks'
if(iprint_DM>2.AND.myid==0) write(io_lun,fmt='(4x,a)') ' Starting Ref To SC Blocks'
! Construct processor ids
if(iprint_DM>1.AND.myid==0) write(io_lun,fmt="(2x,'Scalapack Processor Grid')")
if(iprint_DM>2.AND.myid==0) write(io_lun,fmt="(4x,'Scalapack Processor Grid')")
do ng = 1, proc_groups
n = 1
do i=1,proc_rows
@ -503,7 +503,7 @@ contains
if(n>N_procs_in_pg(ng)) call cq_abort('Ref2SC: Too many processors in group',ng,n)
n = n + 1
end do
if(iprint_DM>1.AND.myid==0) write(io_lun,*) (procid(ng,i,j),j=1,proc_cols)
if(iprint_DM>2.AND.myid==0) write(io_lun,*) (procid(ng,i,j),j=1,proc_cols)
end do
end do
! now build list of blocks and where they go
@ -602,13 +602,13 @@ contains
integer :: i,j,m,n,row,col,proc,ng
integer :: row_max_n, col_max_n, loc_max_row, loc_max_col
if(iprint_DM>1.AND.myid==0) write(io_lun,3)
if(iprint_DM>2.AND.myid==0) write(io_lun,3)
if(iprint_DM>3.AND.myid==0) write(io_lun,4)
! The first row_max_n procs have 1 more block than the rest
row_max_n = mod(blocks_r,proc_rows)
col_max_n = mod(blocks_c,proc_cols)
if(iprint_DM>1.AND.myid==0) write(io_lun,*) 'N for row, col: ',row_max_n, col_max_n
if(iprint_DM>1.AND.myid==0) write(io_lun,*) 'Loc_max_row, col: ',&
if(iprint_DM>2.AND.myid==0) write(io_lun,fmt='(4x,a,2i4)') 'N for row, col: ',row_max_n, col_max_n
if(iprint_DM>2.AND.myid==0) write(io_lun,fmt='(4x,a,2i4)') 'Loc_max_row, col: ',&
aint(real(blocks_r/proc_rows)),aint(real(blocks_c/proc_cols))
my_row = 0
! first record proc_start(:)%rows and %cols map, this is proc proc dependent

116
src/SelfCon_module.f90
View File

@ -178,7 +178,7 @@ contains
use numbers
use global_module, only: iprint_SC, flag_self_consistent, &
flag_SCconverged, IPRINT_TIME_THRES2, &
nspin, spin_factor
nspin, spin_factor, min_layer
use H_matrix_module, only: get_H_matrix
use DMMin, only: FindMinDM
use energy, only: get_energy
@ -215,7 +215,6 @@ contains
call start_backtrace(t=backtrace_timer,who='new_SC_potl',where=area,&
level=backtrace_level,echo=.true.)
!****lat>$
call start_timer(tmr_std_chargescf)
! Build H matrix *with NL and KE*
@ -227,6 +226,7 @@ contains
call start_timer(tmr_l_tmp1,WITH_LEVEL)
if (.not. flag_self_consistent) then
call stop_timer(tmr_std_chargescf)
min_layer = min_layer - 1
call FindMinDM(n_L_iterations, vary_mu, L_tol, &
reset_L, .false.)
call start_timer(tmr_std_chargescf)
@ -234,14 +234,15 @@ contains
call stop_print_timer(tmr_l_tmp1, "new_SC_potl (except H build)", &
IPRINT_TIME_THRES2)
call stop_timer(tmr_std_chargescf)
min_layer = min_layer + 1
return
end if
if (inode == ionode .and. iprint_SC > 0) &
if (inode == ionode .and. iprint_SC + min_layer > 0) &
write (io_lun, &
fmt='(4x,"Starting self-consistency. Tolerance: ",e12.5,/)') &
self_tol
if (record) then
if (inode == ionode .and. iprint_SC > 1) &
if (inode == ionode .and. iprint_SC + min_layer > 1) &
write (io_lun, *) 'Original tol: ', L_tol
SC_tol = self_tol
DMM_tol = L_tol
@ -283,7 +284,7 @@ contains
if (.not. done) then ! Late stage strategy
call start_timer(tmr_l_tmp1, WITH_LEVEL)
reset_L = .true. !.false.
if (inode == ionode .and. iprint_SC > 0) &
if (inode == ionode .and. iprint_SC + min_layer > 0) &
write (io_lun, *) '********** LateSC **********'
call lateSC(record, done, ndone, SC_tol, reset_L, &
@ -306,7 +307,7 @@ contains
call start_timer(tmr_l_tmp1, WITH_LEVEL)
if (record) then ! Fit the C and beta coefficients
if (inode == ionode .and. iprint_SC > 1) then
if (inode == ionode .and. iprint_SC + min_layer > 1) then
write (io_lun, *) ' List of residuals and energies'
if (ndone > 0) then
do i = 1, ndone
@ -315,7 +316,7 @@ contains
end if
end if
call fit_coeff(SCC, SCBeta, SCE, SCR, ndone)
if (inode == ionode .and. iprint_SC > 1) &
if (inode == ionode .and. iprint_SC + min_layer > 1) &
write (io_lun, 6) SCC, SCBeta
end if
@ -913,13 +914,13 @@ contains
use GenBlas
use dimens, only: n_my_grid_points, grid_point_volume
use GenComms, only: gsum, cq_abort, inode, ionode, my_barrier
use io_module, only: dump_charge
use io_module, only: dump_charge, return_prefix
use hartree_module, only: kerker, kerker_and_wdmetric, wdmetric
use global_module, only: ne_in_cell, iprint_SC, area_SC, &
flag_continue_on_SC_fail, &
flag_SCconverged, &
flag_fix_spin_population, nspin, &
spin_factor
spin_factor, min_layer
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl
use maxima_module, only: maxngrid
use store_matrix, only: dump_pos_and_matrices, unit_SCF_save
@ -955,7 +956,10 @@ contains
type(cq_timer) :: backtrace_timer
integer :: backtrace_level
character(len=12) :: subname = "PulayMixSC: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
!****lat<$
if ( present(level) ) backtrace_level = level+1
if ( .not. present(level) ) backtrace_level = -10
@ -971,24 +975,24 @@ contains
R_pul = zero
! write out start information
if (inode == ionode .and. iprint_SC>0) then
write (io_lun, '(8x,a,f6.3,a,f6.3)') &
'Starting Pulay mixing, A_up = ', A(1), ' A_dn = ', A(2)
if (inode == ionode .and. iprint_SC + min_layer>0) then
write (io_lun, '(4x,a,f6.3,a,f6.3)') &
trim(prefix)//' Starting Pulay mixing, A_up = ', A(1), ' A_dn = ', A(2)
if (nspin == 2) then
if (flag_fix_spin_population) then
write (io_lun, '(8x,"Spin populations are fixed.")')
write (io_lun, '(4x,a)') trim(prefix)//" Spin populations are fixed."
else
write (io_lun, '(8x,"Spin populations are to be relaxed.")')
write (io_lun, '(4x,a)') trim(prefix)//" Spin populations are to be relaxed."
end if
else
write (io_lun, '(8x,"Spin non-polarised calculation.")')
write (io_lun, '(4x,a)') trim(prefix)//" Spin non-polarised calculation."
end if
if (flag_Kerker) &
write (io_lun, '(10x,a,f6.3)') &
'with Kerker preconditioning, q0 = ', q0
write (io_lun, '(4x,a,f6.3)') &
trim(prefix)//' with Kerker preconditioning, q0 = ', q0
if (flag_wdmetric) &
write (io_lun, '(10x,a,f6.3)') &
'with wave dependent metric, q1 = ', q1
write (io_lun, '(4x,a,f6.3)') &
trim(prefix)//' with wave dependent metric, q1 = ', q1
end if
! set counters
@ -1042,23 +1046,24 @@ contains
end if
! print residual information
if (inode==ionode) then
if (iprint_SC>1) then
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' RMS residual: ', RA
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' absolute residual (tot) : ', RB
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' absolute residual (frac): ', RC
else if(iprint_SC>0) then
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' residual: ', R0
if (iprint_SC + min_layer>1) then
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' RMS residual: ', RA
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' absolute residual (tot) : ', RB
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' absolute residual (frac): ', RC
else if(iprint_SC + min_layer>0) then
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' residual: ', R0
end if
end if
! check if they have reached tolerance
if (R0 < self_tol .AND. iter >= minitersSC) then ! If we've done minimum number
if (inode == ionode .and. iprint_SC>=0) &
write (io_lun,fmt='(4x,"Reached SCF tolerance of ",e12.5, &
&" after ",i6," iterations")') R0, iter
if (inode == ionode .and. iprint_SC + min_layer>=0) &
write (io_lun,fmt='(4x,a,e12.5, &
&" after ",i6," iterations")') trim(prefix)//" Reached SCF tolerance of ", &
R0, iter
done = .true.
call deallocate_PulayMiXSC_spin
return
@ -1088,7 +1093,9 @@ contains
do iter = 2, maxitersSC
! print out charge
if (flag_DumpChargeDensity .and. iprint_SC > 2) then
! 2019Dec30 tsuyoshi: flag_DumpChargeDensity is introduced, but ...
! is it okay with iprint_SC + min_layer > 1 ?
if (flag_DumpChargeDensity .and. iprint_SC + min_layer > 1) then
if (nspin == 1) then
rho_tot(:) = spin_factor * rho(:,1)
call dump_charge(rho_tot, n_my_grid_points, inode, spin=0)
@ -1144,34 +1151,24 @@ contains
end if
! print residual information
if (inode==ionode) then
if (iprint_SC>1) then
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' RMS residual: ', RA
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' absolute residual (tot) : ', RB
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' absolute residual (frac): ', RC
else if(iprint_SC>0) then
write (io_lun, '(8x,a,i5,a,e12.5)') &
'Pulay iteration ', iter, ' residual: ', R0
if (iprint_SC + min_layer>1) then
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' RMS residual: ', RA
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' absolute residual (tot) : ', RB
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' absolute residual (frac): ', RC
else if(iprint_SC + min_layer>0) then
write (io_lun, '(4x,a,i5,a,e12.5)') &
trim(prefix)//' Pulay iteration ', iter, ' residual: ', R0
end if
end if
! check if they have reached tolerance
if (R0 < self_tol .AND. iter >= minitersSC) then ! Passed minimum number of iterations
! print out charge
if (flag_DumpChargeDensity) then
if (nspin == 1) then
rho_tot(:) = spin_factor * rho(:,1)
call dump_charge(rho_tot, n_my_grid_points, inode, spin=0)
else
call dump_charge(rho(:,1), n_my_grid_points, inode, spin=1)
call dump_charge(rho(:,2), n_my_grid_points, inode, spin=2)
end if
end if
if (inode == ionode .and. iprint_SC>=0) &
write (io_lun,fmt='(4x,"Reached SCF tolerance of ",e12.5, &
&" after ",i6," iterations")') R0, iter
if (inode == ionode .and. iprint_SC + min_layer>=0) &
write (io_lun,fmt='(4x,a,e12.5, &
&" after ",i6," iterations")') trim(prefix)//" Reached SCF tolerance of ", &
R0, iter
done = .true.
call deallocate_PulayMiXSC_spin
return
@ -1187,8 +1184,9 @@ contains
if (R0 > R0_old) &
icounter_fail = icounter_fail + 1
if (icounter_fail > mx_fail) then
if (inode == ionode .and. iprint_SC>0) &
write (io_lun, *) ' Pulay iteration is reset !! at ', iter, &
if (inode == ionode .and. iprint_SC + min_layer>0) &
write (io_lun, fmt='(4x,a,i3,a)') &
trim(prefix)//' Pulay iteration is reset !! at ', iter, &
' th iteration'
reset_Pulay = .true.
end if

23
src/XC_LibXC_v4_module.f90
View File

@ -121,11 +121,7 @@ contains
flag_use_libxc = .true.
call xc_f90_version(vmajor, vminor, vmicro)
if(inode==ionode.AND.iprint_ops>0) then
if(vmajor>2) then
write(io_lun,'("LibXC version: ",I1,".",I1,".",I2)') vmajor, vminor, vmicro
else
write(io_lun,'("LibXC version: ",I1,".",I1)') vmajor, vminor
end if
write(io_lun,'(4x,"LibXC version: ",I2,".",I2,".",I2)') vmajor, vminor, vmicro
end if
! Identify the functional
if(-flag_functional_type<1000) then ! Only exchange OR combined exchange-correlation
@ -195,22 +191,31 @@ contains
if(iprint_ops>2) then
if(vmajor>2) then
write(io_lun,'("The functional ", a, " is ", a, ", it belongs to the ", a, &
write(io_lun,'(4x,"The functional ", a, " is ", a, ", it belongs to the ", a, &
& " family and is defined in the reference(s):")') &
trim(name), trim(kind), trim(family)
j = 0
call xc_f90_info_refs(xc_info(i), j, ref)
do while(j >= 0)
write(io_lun, '(a,i1,2a)') '[', j, '] ', trim(ref)
write(io_lun, '(4x,a,i1,2a)') '[', j, '] ', trim(ref)
call xc_f90_info_refs(xc_info(i), j, ref)
end do
else
write(io_lun,'("The functional ", a, " is ", a, ", and it belongs to the ", a, &
write(io_lun,'(4x,"The functional ", a, " is ", a, ", and it belongs to the ", a, &
& " family")') &
trim(name), trim(kind), trim(family)
end if
else if(iprint_ops>0) then
write(io_lun,'(2x,"Using the ",a," functional ",a)') trim(family),trim(name)
select case(xc_f90_info_kind(xc_info(i)))
case (XC_EXCHANGE)
write(io_lun,'(4x,"Using the ",a," X functional ",a)') trim(family),trim(name)
case (XC_CORRELATION)
write(io_lun,'(4x,"Using the ",a," C functional ",a)') trim(family),trim(name)
case (XC_EXCHANGE_CORRELATION)
write(io_lun,'(4x,"Using the ",a," XC functional ",a)') trim(family),trim(name)
case default
write(io_lun,'(4x,"Using the ",a," functional ",a)') trim(family),trim(name)
end select
else
select case(xc_f90_info_kind(xc_info(i)))
case (XC_EXCHANGE)

24
src/XC_LibXC_v5_module.f90
View File

@ -126,9 +126,9 @@ contains
call xc_f90_version(vmajor, vminor, vmicro)
if(inode==ionode.AND.iprint_ops>0) then
if(vmajor>2) then
write(io_lun,'("LibXC version: ",I1,".",I1,".",I2)') vmajor, vminor, vmicro
write(io_lun,'("LibXC version: ",I2,".",I2,".",I2)') vmajor, vminor, vmicro
else
write(io_lun,'("LibXC version: ",I1,".",I1)') vmajor, vminor
write(io_lun,'("LibXC version: ",I2,".",I2)') vmajor, vminor
end if
end if
! Identify the functional
@ -202,23 +202,11 @@ contains
end select
if(iprint_ops>2) then
if(vmajor>2) then
write(io_lun,'("The functional ", a, " is ", a, ", it belongs to the ", a, &
& " family and is defined in the reference(s):")') &
trim(name), trim(kind), trim(family)
j = 0
ref = xc_f90_func_reference_get_ref(xc_f90_func_info_get_references(xc_info(i),j))
do while(j >= 0)
write(io_lun, '(a,i1,2a)') '[', j, '] ', trim(ref)
ref = xc_f90_func_reference_get_ref(xc_f90_func_info_get_references(xc_info(i),j))
end do
else
write(io_lun,'("The functional ", a, " is ", a, ", and it belongs to the ", a, &
& " family")') &
trim(name), trim(kind), trim(family)
end if
write(io_lun,'(4x,"The functional ", a, " is ", a, ", and it belongs to the ", a, &
" family")') &
trim(name), trim(kind), trim(family)
else if(iprint_ops>0) then
write(io_lun,'(2x,"Using the ",a," functional ",a)') trim(family),trim(name)
write(io_lun,'(4x,"Using the ",a," functional ",a)') trim(family),trim(name)
else
select case(xc_f90_info_kind(xc_info(i)))
case (XC_EXCHANGE)

4
src/block_module.f90
View File

@ -358,8 +358,8 @@ contains
blocks%ng_on_node(nnd) = blocks_per_proc(nnd)
if(blocks_per_proc(nnd)<minblocks) minblocks = blocks_per_proc(nnd)
enddo
if(inode==ionode.AND.iprint_index>0) then
write(io_lun,fmt='(10x,"Minimum blocs/proc is ",i8,". Maximum blocs/proc is ",i8)') minblocks, maxblocks
if(inode==ionode.AND.iprint_index>1) then
write(io_lun,fmt='(4x,"Minimum blocs/proc is ",i8,". Maximum blocs/proc is ",i8)') minblocks, maxblocks
end if
!--- index of first block on current node
blocks%inode_beg=0

427
src/control.f90
View File

@ -140,7 +140,7 @@ contains
use minimise, only: get_E_and_F
use global_module, only: runtype, flag_self_consistent, &
flag_out_wf, flag_write_DOS, wf_self_con, &
flag_opt_cell, optcell_method
flag_opt_cell, optcell_method, min_layer
use input_module, only: leqi
use store_matrix, only: dump_pos_and_matrices
@ -176,23 +176,20 @@ contains
flag_ff, flag_wf, level=backtrace_level)
!
else if ( leqi(runtype, 'cg') ) then
if (flag_opt_cell) then
select case(optcell_method)
case(1)
if (flag_opt_cell) then
select case(optcell_method)
case(1)
call cell_cg_run(fixed_potential, vary_mu, total_energy)
case(2)
case(2)
call full_cg_run_double_loop(fixed_potential, vary_mu, &
total_energy)
case(3)
total_energy)
case(3)
call full_cg_run_single_vector(fixed_potential, vary_mu, &
total_energy)
case(4)
call full_cg_run_double_loop_alt(fixed_potential, vary_mu, &
total_energy)
end select
else
call cg_run(fixed_potential, vary_mu, total_energy)
end if
total_energy)
end select
else
call cg_run(fixed_potential, vary_mu, total_energy)
end if
!
else if ( leqi(runtype, 'md') ) then
call md_run(fixed_potential, vary_mu, total_energy)
@ -279,7 +276,7 @@ contains
y_atom_cell, z_atom_cell, id_glob, &
atom_coord, &
area_general, iprint_MD, &
IPRINT_TIME_THRES1
IPRINT_TIME_THRES1, min_layer
use group_module, only: parts
use minimise, only: get_E_and_F
use move_atoms, only: adapt_backtrack_linemin, safemin2, cg_line_min, safe, backtrack
@ -335,6 +332,7 @@ contains
energy1 = zero
dE = zero
! Find energy and forces
min_layer = min_layer - 1
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., .true.)
call dump_pos_and_matrices
call get_maxf(max)
@ -342,6 +340,7 @@ contains
write(io_lun,'(/4x,"GeomOpt - Iter: ",i4," MaxF: ",f12.8," E: ",e16.8," dE: ",f12.8/)') &
0, max, energy0, dE
end if
min_layer = min_layer + 1
iter = 1
ggold = zero
@ -412,6 +411,7 @@ contains
end if
old_force = tot_force
! Minimise in this direction
!min_layer = min_layer - 1
if(cg_line_min==safe) then
call safemin2(x_new_pos, y_new_pos, z_new_pos, cg, energy0,&
energy1, fixed_potential, vary_mu, energy1)
@ -427,6 +427,7 @@ contains
! Analyse forces
g0 = dot(length, tot_force, 1, tot_force, 1)
call get_maxf(max)
!min_layer = min_layer + 1
! Output and energy changes
dE = energy1 - energy0
@ -593,7 +594,7 @@ contains
flag_move_atom,rcellx, rcelly, rcellz, &
flag_Multisite,flag_SFcoeffReuse, &
atom_coord, flag_quench_MD, atomic_stress, &
non_atomic_stress, flag_heat_flux
non_atomic_stress, flag_heat_flux, min_layer
use group_module, only: parts
use minimise, only: get_E_and_F
use move_atoms, only: velocityVerlet, updateIndices, &
@ -704,11 +705,13 @@ contains
call mdl%get_cons_qty
! Find energy and forces
min_layer = min_layer - 1
if (flag_fire_qMD) then
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., .true.)
else
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., .false.)
end if
min_layer = min_layer + 1
mdl%dft_total_energy = energy0
! XL-BOMD
@ -805,6 +808,7 @@ contains
! may change after the atomic positions are updated.
call check_move_atoms(flag_movable)
min_layer = min_layer - 1
if (flag_fire_qMD) then
call get_E_and_F(fixed_potential, vary_mu, energy1, .true., .true.,iter)
call check_stop(done, iter) !2019/Nov/14
@ -822,6 +826,7 @@ contains
call dump_pos_and_matrices(index=0,MDstep=iter,velocity=ion_velocity)
call vVerlet_v_dthalf(MDtimestep,ion_velocity,tot_force,flag_movable,second_call)
end if
min_layer = min_layer + 1
! Update DFT energy
mdl%dft_total_energy = energy1
!******
@ -1935,7 +1940,7 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
y_atom_cell, z_atom_cell, id_glob, &
atom_coord, area_general, flag_pulay_simpleStep, &
flag_diagonalisation, nspin, flag_LmatrixReuse, &
flag_SFcoeffReuse
flag_SFcoeffReuse, min_layer
use group_module, only: parts
use minimise, only: get_E_and_F
use move_atoms, only: pulayStep, velocityVerlet, &
@ -2010,10 +2015,12 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
energy1 = zero
dE = zero
! Find energy and forces
min_layer = min_layer - 1
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., &
.false.)
call dump_pos_and_matrices
call get_maxf(max)
min_layer = min_layer + 1
iter = 0
ggold = zero
energy1 = energy0
@ -2558,7 +2565,7 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
atom_coord, rcellx, rcelly, rcellz, &
area_general, iprint_MD, &
IPRINT_TIME_THRES1, cell_en_tol, &
cell_constraint_flag, cell_stress_tol
cell_constraint_flag, cell_stress_tol, min_layer
use group_module, only: parts
use minimise, only: get_E_and_F
use move_atoms, only: safemin_cell, enthalpy, enthalpy_tolerance
@ -2609,7 +2616,9 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
energy1 = zero
dE = zero
! Find energy and forces
min_layer = min_layer - 1
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., .true.)
min_layer = min_layer + 1
iter = 1
reset_iter = 1
ggold = zero
@ -2687,9 +2696,11 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
! Minimise in this direction. Constraint information is also used within
! safemin_cell. Look in move_atoms.module.f90 for further information.
min_layer = min_layer - 1
call safemin_cell(new_rcellx, new_rcelly, new_rcellz, search_dir_x, &
search_dir_y, search_dir_z, search_dir_mean, press, &
enthalpy0, enthalpy1, fixed_potential, vary_mu)
min_layer = min_layer + 1
! Output positions to UpdatedAtoms.dat
if (myid == 0 .and. iprint_gen > 1) then
do i = 1, ni_in_cell
@ -3135,377 +3146,6 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
end subroutine full_cg_run_double_loop
!!***
! Keeping this temporarily so that we can test the relative efficiency of
! alternating full ionic and full cell optimisation (full_cg_run_double_loop)
! and full ionic with single line minimisation cell optimisation (this routine)
! Use cell optimisation method 4 for this
subroutine full_cg_run_double_loop_alt(fixed_potential, vary_mu, total_energy)
! Module usage
use numbers
use units
use global_module, only: iprint_gen, ni_in_cell, x_atom_cell, &
y_atom_cell, z_atom_cell, id_glob, &
atom_coord, area_general, iprint_MD, &
IPRINT_TIME_THRES1, &
cell_en_tol, cell_stress_tol, &
rcellx, rcelly, rcellz
use group_module, only: parts
use minimise, only: get_E_and_F
use move_atoms, only: safemin2, safemin_cell, enthalpy, &
enthalpy_tolerance
use GenComms, only: inode, ionode
use GenBlas, only: dot
use force_module, only: tot_force, stress
use io_module, only: write_atomic_positions, pdb_template, &
check_stop, write_xsf
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl
use timer_module
use store_matrix, only: dump_InfoMatGlobal, dump_pos_and_matrices
use md_control, only: flag_write_xsf, target_pressure
implicit none
! Passed variables
! Shared variables needed by get_E_and_F for now (!)
logical :: vary_mu, fixed_potential
real(double) :: total_energy
! Local variables for ionic relaxation
real(double) :: energy0, energy1, max, g0, dE, gg, ggold, gamma,gg1, &
enthalpy0, enthalpy1, dH, press
integer :: i,j,k,iter,length, jj, lun, stat
logical :: done_ions, done_cell
type(cq_timer) :: tmr_l_iter
real(double), allocatable, dimension(:,:) :: cg, old_force
real(double), allocatable, dimension(:) :: x_new_pos, y_new_pos,&
z_new_pos
! Local variables for cell relaxation
real(double) :: new_rcellx, new_rcelly, new_rcellz, search_dir_x, &
search_dir_y, search_dir_z, stressx, stressy, stressz, &
RMSstress, newRMSstress, dRMSstress, search_dir_mean, &
mean_stress, energy2, volume, stress_diff, max_stress, stress_target
integer :: iter_cell
if (inode==ionode .and. iprint_MD > 2) &
write(io_lun,'(2x,a)') "control/full_cg_run_double_loop"
allocate(cg(3,ni_in_cell), STAT=stat)
if (stat /= 0) &
call cq_abort("Error allocating cg in control: ",&
ni_in_cell, stat)
allocate(old_force(3,ni_in_cell), STAT=stat)
allocate(x_new_pos(ni_in_cell), y_new_pos(ni_in_cell), &
z_new_pos(ni_in_cell), STAT=stat)
if (stat/=0) &
call cq_abort("Error allocating _new_pos in control: ", &
ni_in_cell,stat)
call reg_alloc_mem(area_general, 6 * ni_in_cell, type_dbl)
search_dir_z = zero
search_dir_x = zero
search_dir_y = zero
search_dir_mean = zero
cg = zero
done_cell = .false.
length = 3 * ni_in_cell
if (inode==ionode .and. iprint_gen > 0) then
write(io_lun,'(4x,"Welcome to cg_run. Doing ",i4," steps")') MDn_steps
write(io_lun,'(4x,"Force tolerance: ",f20.10)') MDcgtol
write(io_lun,'(4x,"Enthalpy tolerance: ",f20.10)') enthalpy_tolerance
end if
energy0 = total_energy
energy1 = zero
energy2 = zero
dE = zero
! Find energy and forces
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., .true.)
call dump_pos_and_matrices
call get_maxf(max)
press = target_pressure/HaBohr3ToGPa
! Stress tolerance in Ha/Bohr3
stress_target = cell_stress_tol/HaBohr3ToGPa
enthalpy0 = enthalpy(energy0, press)
dH = zero
if (inode==ionode) then
write(io_lun,'(/4x,"GeomOpt - Iter: ",i4," MaxF: ",f12.8," H: ", e16.8," dH: ",f12.8/)') &
0, max, enthalpy0, dH
end if
if (inode == ionode .and. iprint_gen > 0) then
write (io_lun, fmt='(/4x,"Starting full cell optimisation"/)')
write(io_lun,*) "Initial cell dims", rcellx, rcelly, rcellz
end if
iter = 1
iter_cell = 1
! Cell loop
do while (.not. done_cell)
<<<<<<< HEAD
call cg_run(fixed_potential, vary_mu, energy1)
energy0 = energy1
=======
ggold = zero
old_force = zero
energy1 = energy0
enthalpy1 = enthalpy0
done_ions = .false.
do while (.not. done_ions) ! ionic loop
call start_timer(tmr_l_iter, WITH_LEVEL)
! Construct ratio for conjugacy
gg = zero
gg1 = zero! PR
do j = 1, ni_in_cell
gg = gg + &
tot_force(1,j) * tot_force(1,j) + &
tot_force(2,j) * tot_force(2,j) + &
tot_force(3,j) * tot_force(3,j)
gg1 = gg1 + &
tot_force(1,j) * old_force(1,j) + &
tot_force(2,j) * old_force(2,j) + &
tot_force(3,j) * old_force(3,j)
end do
if (abs(ggold) < 1.0e-6_double) then
gamma = zero
else
gamma = (gg-gg1)/ggold ! PR - change to gg/ggold for FR
end if
if(gamma<zero) gamma = zero
if (inode == ionode .and. iprint_MD > 2) &
write (io_lun,*) ' CHECK :: Force Residual = ', &
for_conv * sqrt(gg)/ni_in_cell
if (inode == ionode .and. iprint_MD > 2) &
write (io_lun,*) ' CHECK :: gamma = ', gamma
if (CGreset) then
if (gamma > one) then
if (inode == ionode) &
write(io_lun,*) ' CG direction is reset! '
gamma = zero
end if
end if
if (inode == ionode) &
write (io_lun, fmt='(/4x,"Atomic relaxation CG iteration: ",i5)') iter
ggold = gg
! Build search direction
do j = 1, ni_in_cell
jj = id_glob(j)
cg(1,j) = gamma*cg(1,j) + tot_force(1,jj)
cg(2,j) = gamma*cg(2,j) + tot_force(2,jj)
cg(3,j) = gamma*cg(3,j) + tot_force(3,jj)
x_new_pos(j) = x_atom_cell(j)
y_new_pos(j) = y_atom_cell(j)
z_new_pos(j) = z_atom_cell(j)
end do
old_force = tot_force
! Minimise in this direction
call safemin2(x_new_pos, y_new_pos, z_new_pos, cg, energy0, &
energy1, fixed_potential, vary_mu, energy1)
! Output positions
if (inode==ionode .and. iprint_gen > 1) then
write(io_lun,'(4x,a4,a15)') "Atom", "Position"
do i = 1, ni_in_cell
write(io_lun,'(4x,i8,3f15.8)') i,atom_coord(:,i)
end do
end if
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
if (flag_write_xsf) call write_xsf('trajectory.xsf', iter)
! Analyse forces
g0 = dot(length, tot_force, 1, tot_force, 1)
call get_maxf(max)
! Output and energy changes
enthalpy0 = enthalpy(energy0, press)
enthalpy1 = enthalpy(energy1, press)
dE = energy1 - energy0
dH = enthalpy1 - enthalpy0
energy0 = energy1
if (inode==ionode) then
write(io_lun,'(/4x,"GeomOpt - Iter: ",i4," MaxF: ",f12.8," H: ",e16.8," dH: ",f12.8/)') &
iter, max, enthalpy1, en_conv*dH
if (iprint_MD > 1) then
write(io_lun,'(4x,"Force Residual: ",f20.10," ",a2,"/",a2)') &
for_conv*sqrt(g0/ni_in_cell), en_units(energy_units), &
d_units(dist_units)
write(io_lun,'(4x,"Maximum force: ",f20.10)') max
write(io_lun,'(4x,"Force tolerance: ",f20.10)') MDcgtol
write(io_lun,'(4x,"Maximum stress ",3f10.6)') &
max_stress
write(io_lun,'(4x,"Stress tolerance: ",f20.10)') &
cell_stress_tol
write(io_lun,'(4x,"Enthalpy change: ",f20.10," ",a2)') &
en_conv*dH, en_units(energy_units)
write(io_lun,'(4x,"Enthalpy tolerance: ",f20.10)') &
enthalpy_tolerance
end if
end if
iter = iter + 1
if (iter > MDn_steps) then
done_ions = .true.
if (inode==ionode) &
write (io_lun, fmt='(4x,"Exceeded number of CG steps: ",i6)') &
iter
end if
if (abs(max) < MDcgtol) then
done_ions = .true.
if (inode==ionode) &
write (io_lun, &
fmt='(4x,"Maximum force below threshold: ",f12.5)') max
end if
call dump_pos_and_matrices
call stop_print_timer(tmr_l_iter, "a CG iteration", IPRINT_TIME_THRES1)
if (.not. done_ions) call check_stop(done_ions, iter)
end do ! ionic loop
>>>>>>> Tidying output for CG searches
enthalpy0 = enthalpy(energy0, press)
volume = rcellx*rcelly*rcellz
stressx = -stress(1,1)!/volume
stressy = -stress(2,2)!/volume
stressz = -stress(3,3)!/volume
mean_stress = (stressx + stressy + stressz)/3
RMSstress = sqrt(((stressx*stressx) + (stressy*stressy) + (stressz*stressz))/3)
gamma = zero ! steepest descent only for cell iteration
if (inode == ionode .and. iprint_gen > 0) &
write (io_lun, fmt='(/4x,"Lattice vector relaxation iteration: ",i5)') iter_cell
ggold = gg
!Build search direction.
! If the volume constraint is set, there is only one search direction!
! This is the direction which minimises the mean stress.
! if (leqi(cell_constraint_flag, 'volume')) then
! search_dir_mean = gamma*search_dir_mean + mean_stress
! else
search_dir_x = gamma*search_dir_x + stressx - press*volume
search_dir_y = gamma*search_dir_y + stressy - press*volume
search_dir_z = gamma*search_dir_z + stressz - press*volume
! end if
new_rcellx = rcellx
new_rcelly = rcelly
new_rcellz = rcellz
! Minimise in this direction. Constraint information is also used within
! safemin_cell. Look in move_atoms.module.f90 for further information.
call safemin_cell(new_rcellx, new_rcelly, new_rcellz, search_dir_x, &
search_dir_y, search_dir_z, search_dir_mean, press, &
enthalpy0, enthalpy1, fixed_potential, vary_mu)
! Output positions to UpdatedAtoms.dat
if (inode==ionode .and. iprint_MD > 1) then
write(io_lun,'(4x,a4,a15)') "Atom", "Position"
do i = 1, ni_in_cell
write(io_lun,'(4x,i8,3f15.8)') i,atom_coord(:,i)
end do
end if
call write_atomic_positions("UpdatedAtoms.dat", trim(pdb_template))
! Analyse Stresses and energies
dH = enthalpy0 - enthalpy1
call get_maxf(max)
newRMSstress = sqrt(((stress(1,1)*stress(1,1)) + &
(stress(2,2)*stress(2,2)) + &
(stress(3,3)*stress(3,3)))/three)
dRMSstress = RMSstress - newRMSstress
volume = rcellx*rcelly*rcellz
max_stress = zero
do i=1,3
stress_diff = abs(press*volume + stress(i,i))/volume
if (stress_diff > max_stress) max_stress = stress_diff
end do
if (inode==ionode) then
write(io_lun,'(/4x,"GeomOpt + Iter: ",i4," MaxF: ",f12.8," H: " ,e16.8," dH: ",f12.8/)') &
iter, max, enthalpy1, en_conv*dH
if (iprint_MD > 1) then
write(io_lun,'(4x,"Force Residual: ",f20.10," ",a2,"/",a2)') &
for_conv*sqrt(g0/ni_in_cell), en_units(energy_units), &
d_units(dist_units)
write(io_lun,'(4x,"Maximum force: ",f20.10)') max
write(io_lun,'(4x,"Force tolerance: ",f20.10)') MDcgtol
write(io_lun,'(4x,"Maximum stress ",e14.6," Ha/Bohr**3")') &
max_stress
write(io_lun,'(4x,"Simulation cell volume ",e14.6," Bohr**3")') volume
write(io_lun,'(4x,"Maximum stress ",f14.6," GPa")') &
max_stress*HaBohr3ToGPa
write(io_lun,'(4x,"Stress tolerance: ",f14.6," GPa")') &
cell_stress_tol
write(io_lun,'(4x,"Enthalpy change: ",f20.10," ",a2)') &
en_conv*dH, en_units(energy_units)
write(io_lun,'(4x,"Enthalpy tolerance: ",f20.10)') &
enthalpy_tolerance
end if
end if
iter = iter + 1
iter_cell = iter_cell + 1
energy0 = energy1
enthalpy0 = enthalpy1
! Check exit criteria
if (iter > MDn_steps) then
done_cell = .true.
if (inode==ionode) &
write (io_lun, fmt='(4x,"Exceeded number of SD steps: ",i4)') iter
end if
! If the force convergence criterion has been met after this step,
! we can assume the ionic positions are correct, THEN check cell
! convergence
if (abs(dH) < enthalpy_tolerance) then
if (abs(max) < MDcgtol .and. max_stress < stress_target) then
done_cell = .true.
if (inode==ionode) then
write(io_lun,'(/4x,a,i4,a,i4,a)') "GeomOpt converged in ", &
iter-iter_cell," ionic steps and ", iter_cell, " cell steps"
write(io_lun, fmt='(4x,"Maximum force: ",f20.10)') max
write(io_lun, fmt='(4x,"Force tolerance: ",f20.10)') MDcgtol
write(io_lun, fmt='(4x,"Enthalpy change: ",f20.10)') dH
write(io_lun, fmt='(4x,"Enthalpy tolerance: ",f20.10)') &
enthalpy_tolerance
write(io_lun, fmt='(4x,"Maximum stress: ",f20.10," GPa")') &
max_stress*HaBohr3ToGPa
write(io_lun, fmt='(4x,"Stress tolerance: ",f20.10," GPa")') &
cell_stress_tol
end if
end if
else
continue
end if
call stop_print_timer(tmr_l_iter, "a CG iteration", IPRINT_TIME_THRES1)
if (.not. done_cell) call check_stop(done_cell, iter_cell)
end do
if (inode == ionode .and. iprint_gen > 0) then
write(io_lun, fmt='(/4x,"Finished full cell optimisation"/)')
write(io_lun,*) "Final cell dims", rcellx, rcelly, rcellz
end if
deallocate(z_new_pos, y_new_pos, x_new_pos, STAT=stat)
if (stat /= 0) &
call cq_abort("Error deallocating _new_pos in control: ", &
ni_in_cell,stat)
deallocate(old_force, STAT=stat)
deallocate(cg, STAT=stat)
if (stat /= 0) &
call cq_abort("Error deallocating cg in control: ", ni_in_cell,stat)
call reg_dealloc_mem(area_general, 6*ni_in_cell, type_dbl)
3 format(4x,'*** CG step ',i4,' Gamma: ',f14.8)
4 format(4x,'Enthalpy change: ',f15.8,' ',a2)
5 format(4x,'Force Residual: ',f15.10,' ',a2,'/',a2)
6 format(4x,'Maximum force component: ',f15.8,' ',a2,'/',a2)
7 format(4x,3f15.8)
end subroutine full_cg_run_double_loop_alt
!!****f* control/full_cg_run_single_vector *
!!
!! NAME
@ -3542,7 +3182,7 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
atom_coord, rcellx, rcelly, rcellz, &
area_general, iprint_MD, &
IPRINT_TIME_THRES1, &
cell_stress_tol
cell_stress_tol, min_layer
use group_module, only: parts
use minimise, only: get_E_and_F
use move_atoms, only: safemin_full, cq_to_vector, enthalpy, &
@ -3608,9 +3248,11 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
cell_ref(3) = rcellz
! Find energy and forces
min_layer = min_layer - 1
call get_E_and_F(fixed_potential, vary_mu, energy0, .true., .true.)
call dump_pos_and_matrices
call get_maxf(max)
min_layer = min_layer + 1
enthalpy0 = enthalpy(energy0, press)
if (inode==ionode) then
write(io_lun,'(/4x,"GeomOpt - Iter: ",i4," MaxF: ",f12.8," H: ",e16.8," dH: ",f12.8/)') &
@ -3666,9 +3308,10 @@ subroutine update_pos_and_box(baro, nequil, flag_movable)
end do
force_old = force
! Minimise in this direction
min_layer = min_layer - 1
call safemin_full(config, cg, cell_ref, enthalpy0, enthalpy1, &
press, fixed_potential, vary_mu, enthalpy1)
min_layer = min_layer + 1
! Output positions
if (inode==ionode .and. iprint_gen > 1) then
write(io_lun,'(4x,a4,a15)') "Atom", "Position"

2
src/cover_module.f90
View File

@ -354,7 +354,7 @@ contains
!endif ! members.AND.nm_group>0
enddo
if(members) then ! Now generate member information
if(inode==ionode.AND.iprint_index>1) write(io_lun,*) 'Members in covering set: ',nm_in_cover
if(inode==ionode.AND.iprint_index>2) write(io_lun,fmt='(4x,a,i7)') 'Members in covering set: ',nm_in_cover
if(members) set%icover_ibeg(1)=1
set%mx_mcover = nm_in_cover
call start_timer(tmr_std_allocation)

62
src/density_module.f90
View File

@ -209,7 +209,7 @@ contains
IPRINT_TIME_THRES3, nspin, &
spin_factor, &
flag_fix_spin_population, &
flag_neutral_atom, area_SC
flag_neutral_atom, area_SC, min_layer
use block_module, only: nx_in_block, ny_in_block, &
nz_in_block, n_pts_in_block
use group_module, only: blocks, parts
@ -224,7 +224,8 @@ contains
use maxima_module, only: maxngrid
use species_module, only: charge, charge_up, charge_dn, type_species
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, &
type_dbl
type_dbl
use io_module, only: return_prefix
implicit none
@ -249,6 +250,8 @@ contains
type(cq_timer) :: tmr_l_tmp1
type(cq_timer) :: backtrace_timer
integer :: backtrace_level, stat
character(len=15) :: subname = "set_atom_dens: "
character(len=120) :: prefix
!****lat<$
if ( present(level) ) backtrace_level = level+1
@ -256,10 +259,11 @@ contains
call start_backtrace(t=backtrace_timer,who='set_density', &
where=area,level=backtrace_level)
!****lat>$
prefix = return_prefix(subname, min_layer)
if (inode == ionode .and. iprint_SC >= 2) then
write (io_lun, fmt='(2x,"Entering set_density")')
if (flag_InitialAtomicSpin) write (io_lun, fmt='(2x,"Initial atomic spins are read from input file")')
write (io_lun, fmt='(4x,a)') trim(prefix)//"Entering"
if (flag_InitialAtomicSpin) write (io_lun, fmt='(4x,a)') &
trim(prefix)//"Initial atomic spins are read from input file"
endif
call start_timer(tmr_std_chargescf)
@ -413,6 +417,8 @@ contains
call gsum(grid_electrons)
! Scaling factor for renormalizing atomic density
scale = ne_in_cell/grid_electrons
! Renormalize the atomic density
density_atom(:) = density_atom(:) * scale
!write(*,*) "SCALE in atomic_density = ", scale, ne_in_cell, grid_electrons
@ -431,8 +437,9 @@ contains
density(:,spin) = density(:,spin) * density_scale(spin)
if (inode == ionode .and. iprint_SC > 0) &
write (io_lun, &
fmt='(10x,"In set_atomic_density, electrons (spin=",i1,"): ",f20.12)') &
spin, density_scale(spin) * grid_electrons
fmt='(4x,a,i1,a,f20.12)') &
trim(prefix)//"In set_atomic_density, electrons (spin=", spin, "): ", &
density_scale(spin) * grid_electrons
enddo
else
! -- Scale for spin
@ -441,15 +448,20 @@ contains
! If they are unequal, then scale half*density_atom by relative population
density_scale(spin) = ne_spin_in_cell(spin) / (half * grid_electrons) ! Calculate relative population for spin channel
density(:,spin) = density_scale(spin) * half * density_atom(:) ! Assign appropriate amount of (half*density_atom) to spin channel
if (inode == ionode .and. iprint_SC > 0) &
write (io_lun, &
fmt='(10x,"In set_atomic_density, electrons (spin=",i1,"): ",f20.12)') &
spin, density_scale(spin) * half * grid_electrons
if (inode == ionode .and. iprint_SC > 0) then
if(nspin==1) then
write (io_lun, fmt='(4x,a,f20.12)') &
trim(prefix)//"In set_atomic_density, electrons : ", &
spin_factor * density_scale(spin) * half * grid_electrons
else
write (io_lun, fmt='(4x,a,i1,a,f20.12)') &
trim(prefix)//"In set_atomic_density, electrons (spin=", &
spin, "): ", density_scale(spin) * half * grid_electrons
end if
end if
end do
endif
end if ! flag_set_density
! Renormalize the atomic density
density_atom(:) = density_atom(:) * scale
if(.NOT.flag_neutral_atom) then
deallocate(density_atom, STAT=stat)
@ -506,7 +518,7 @@ contains
use global_module, only: rcellx, rcelly, rcellz, id_glob, &
ni_in_cell, iprint_SC, &
species_glob, dens, ne_in_cell, &
IPRINT_TIME_THRES3
IPRINT_TIME_THRES3, min_layer
use block_module, only: nx_in_block, ny_in_block, &
nz_in_block, n_pts_in_block
use group_module, only: blocks, parts
@ -518,6 +530,7 @@ contains
use dimens, only: n_my_grid_points, grid_point_volume
use GenBlas, only: rsum, scal
use timer_module
use io_module, only: return_prefix
implicit none
@ -537,9 +550,12 @@ contains
real(double) :: pcc_density !P.C.C. charge density interpolated
real(double) :: a, b, c, d, r1, r2, r3, r4, rr
type(cq_timer) :: tmr_l_tmp1
character(len=9) :: subname = "Set PCC: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
if (inode == ionode .and. iprint_SC >= 2) &
write (io_lun, fmt='(2x,"Entering set_density_pcc")')
write (io_lun, fmt='(4x,a)') trim(prefix)//"Entering"
call start_timer(tmr_std_chargescf)
call start_timer(tmr_l_tmp1, WITH_LEVEL)
@ -767,8 +783,9 @@ contains
use GenComms, only: gsum
use global_module, only: iprint_SC, ni_in_cell, &
flag_Becke_weights, nspin, &
spin_factor, sf, paof, atomf
spin_factor, sf, paof, atomf, min_layer
use functions_on_grid, only: gridfunctions, fn_on_grid
use io_module, only: return_prefix
implicit none
@ -783,7 +800,10 @@ contains
type(cq_timer) :: backtrace_timer
integer :: backtrace_level
integer :: blk, i_count_alpha, n, n_i, n_point, spin
character(len=12) :: subname = "get_density: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
!****lat<$
if ( present(level) ) backtrace_level = level+1
if ( .not. present(level) ) backtrace_level = -10
@ -793,8 +813,8 @@ contains
call start_timer(tmr_std_chargescf)
if (inode == ionode .and. iprint_SC >= 2) &
write (io_lun,fmt='(2x,"Entering get_electronic_density")')
if (inode == ionode .and. iprint_SC + min_layer >= 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Entering'
do spin = 1, nspin
! matK->matKatomf backtransformation for contracted SFs
@ -843,9 +863,9 @@ contains
rsum(n_my_grid_points, denout(:,spin), 1)
call gsum(electrons(spin))
if (inode == ionode .and. iprint_SC > 1) &
write (io_lun, '(2x,"Electrons (spin=",i1,"): ",f25.15)') &
spin, electrons(spin)
if (inode == ionode .and. iprint_SC + min_layer > 1) &
write (io_lun, '(4x,a,i1,a,f25.15)') &
trim(prefix)//"Electrons (spin=",spin, "): ", electrons(spin)
end do ! spin
if (flag_Becke_weights) &

147
src/energy_module.f90
View File

@ -161,7 +161,7 @@ contains
flag_perform_cdft, &
flag_vdWDFT, &
flag_exx, exx_alpha, &
flag_neutral_atom
flag_neutral_atom, min_layer
use pseudopotential_common, only: core_correction, &
flag_neutral_atom_projector
use density_module, only: electron_number
@ -257,7 +257,7 @@ contains
if (inode == ionode) then
if(print_Harris) then
!if(iprint_gen>=1) write(io_lun,2) electrons
if (iprint_gen >= 1) then
if (iprint_gen + min_layer >= 1) then
if(flag_neutral_atom) then
write (io_lun, 1) en_conv*band_energy, en_units(energy_units)
write (io_lun,33) en_conv*hartree_energy_drho, en_units(energy_units)
@ -295,18 +295,18 @@ contains
end if
if (abs(entropy) >= RD_ERR) then
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write(io_lun,10) en_conv*total_energy, en_units(energy_units)
if (flag_check_Diag) then
select case (SmearingType)
case (0) ! Fermi smearing
if (entropy < zero) &
call cq_warn(sub_name,'Calculated entropy is less than zero; something is wrong ', entropy)
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun,14) en_conv*(total_energy-half*entropy), &
en_units(energy_units)
case (1) ! Methfessel-Paxton smearing
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun,16) &
en_conv * (total_energy - &
(real(MPOrder+1,double) / &
@ -314,17 +314,17 @@ contains
en_units(energy_units)
end select
else
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun,14) en_conv*(total_energy-half*entropy), &
en_units(energy_units)
end if
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun,15) en_conv*(total_energy-entropy), &
en_units(energy_units)
else
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun,10) en_conv*total_energy, en_units(energy_units)
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun, '(10x,"(TS=0 as O(N) or entropic &
&contribution is negligible)")')
end if
@ -353,16 +353,16 @@ contains
if (flag_perform_cdft) total_energy2 = total_energy2 + cdft_energy
if (flag_dft_d2) total_energy2 = total_energy2 + disp_energy
if (inode == ionode .and. iprint_gen>=1) then
if (inode == ionode .and. iprint_gen + min_layer>=1) then
write(io_lun,13) en_conv*total_energy2, en_units(energy_units)
end if
end if
! print electron number and spin polarisation information
if (iprint_gen >= 0) then
if (iprint_gen + min_layer >= 0) then
call electron_number(electrons)
if (inode == ionode) then
if (iprint_gen >= 1) then
if (iprint_gen + min_layer >= 1) then
electrons_tot = electrons(1) + electrons(nspin)
write (io_lun,18) electrons_tot
if (nspin == 2) then
@ -462,7 +462,9 @@ contains
flag_self_consistent, &
flag_perform_cdft, &
flag_vdWDFT, &
flag_exx, exx_alpha, flag_neutral_atom
flag_exx, exx_alpha, &
flag_neutral_atom, min_layer
use DFT_D2, only: disp_energy
use density_module, only: electron_number
use pseudopotential_common, only: core_correction, &
flag_neutral_atom_projector
@ -565,13 +567,27 @@ contains
if (inode == ionode) then
electrons_tot = electrons(1) + electrons(nspin)
!
if(iprint_gen==0) then
if(iprint_gen + min_layer>=-1) then
if (nspin == 1) then
write (io_lun,fmt='(4x," | Number of electrons = ",f16.6)') electrons_tot
write (io_lun,fmt='(4x,"| Number of electrons = ",f16.6)') electrons_tot
else if(nspin == 2) then
write (io_lun,fmt='(4x," | Number of electrons (u/d)= ",2f16.6)') electrons(1),electrons(nspin)
write (io_lun,fmt='(4x,"| Number of electrons (u/d)= ",2f16.6)') electrons(1),electrons(nspin)
end if
else if (iprint_gen >= 1) then
else if (iprint_gen + min_layer ==1) then
write (io_lun, 6) en_conv * band_energy, en_units(energy_units)
if(flag_neutral_atom) then
write (io_lun,68) en_conv * screened_ion_interaction_energy, en_units(energy_units)
else
write (io_lun, 8) en_conv * ion_interaction_energy, en_units(energy_units)
end if
if(flag_neutral_atom) then
write (io_lun,11) en_conv*(- hartree_energy_drho - hartree_energy_drho_atom_rho), &
en_units(energy_units)
else
write (io_lun,11) en_conv* delta_E_hartree, en_units(energy_units)
end if
write (io_lun,12) en_conv* delta_E_xc, en_units(energy_units)
else if (iprint_gen + min_layer >= 2) then
write (io_lun, *)
!write (io_lun, *)
!write (io_lun, 1)
@ -635,13 +651,13 @@ contains
if (entropy < zero) &
call cq_warn(sub_name, 'Calculated entropy is less than zero; something is wrong ', entropy)
!
if (iprint_gen >= 0) &
if (iprint_gen + min_layer >= 0) &
write (io_lun,14) en_conv*(total_energy1-half*entropy), &
en_units(energy_units)
!
!
case (1) ! Methfessel-Paxton smearing
if (iprint_gen >= 0) &
if (iprint_gen + min_layer >= 0) &
write (io_lun,16) en_conv * (total_energy1 - &
(real(MPOrder+1,double) / &
real(MPOrder+2,double))*entropy), &
@ -651,13 +667,13 @@ contains
end select
!
else
if (iprint_gen >= 0) &
if (iprint_gen + min_layer >= 0) &
write (io_lun,14) en_conv*(total_energy1-half*entropy), &
en_units(energy_units)
end if
!
!
if (iprint_gen >= 1) &
if (iprint_gen + min_layer >= 1) &
write (io_lun,15) en_conv*(total_energy1-entropy), &
en_units(energy_units)
else
@ -702,11 +718,12 @@ contains
hartree_energy_total_rho
if (inode == ionode) then
if (iprint_gen >= 0) then
if (iprint_gen + min_layer >= -1) then
write(io_lun,10) en_conv*total_energy1, en_units(energy_units)
if(iprint_gen>0) then
write(io_lun,13) en_conv*total_energy2, en_units(energy_units)
write(io_lun,22) en_conv*(total_energy1 - total_energy2), &
if(iprint_gen + min_layer>0) then
write(io_lun,13) en_conv*total_energy2, en_units(energy_units)
if(iprint_gen + min_layer>1) &
write(io_lun,22) en_conv*(total_energy1 - total_energy2), &
en_units(energy_units)
end if
end if
@ -723,7 +740,9 @@ contains
if (inode == ionode) then
electrons_tot = electrons(1) + electrons(nspin)
if (iprint_gen >= 1) then
if(iprint_gen + min_layer==1) then
write (io_lun,24) electrons_tot
elseif (iprint_gen + min_layer >= 2) then
write (io_lun,23)
write (io_lun,24) electrons_tot
write (io_lun,25) electrons_tot2
@ -756,51 +775,51 @@ contains
2 format(4x, ' ')
6 format(4x,' |* band energy as 2Tr[K.H] = ',f25.15,' ',a2)
7 format(4x,' | hartree energy (rho) = ',f25.15,' ',a2)
8 format(4x,' | ion-ion energy = ',f25.15,' ',a2)
67 format(4x,' | hartree energy (drho) = ',f25.15,' ',a2)
68 format(4x,' | screened ion-ion energy = ',f25.15,' ',a2)
9 format(4x,' | kinetic energy = ',f25.15,' ',a2)
6 format(4x,'|* band energy as 2Tr[K.H] = ',f25.15,' ',a2)
7 format(4x,'| hartree energy (rho) = ',f25.15,' ',a2)
8 format(4x,'| ion-ion energy = ',f25.15,' ',a2)
67 format(4x,'| hartree energy (drho) = ',f25.15,' ',a2)
68 format(4x,'| screened ion-ion energy = ',f25.15,' ',a2)
9 format(4x,'| kinetic energy = ',f25.15,' ',a2)
30 format(4x,' |* xc total energy = ',f25.15,' ',a2)
31 format(4x,' | DFT exchange = ',f25.15,' ',a2)
32 format(4x,' | DFT correlation = ',f25.15,' ',a2)
33 format(4x,' | EXX contribution = ',f25.15,' ',a2)
30 format(4x,'|* xc total energy = ',f25.15,' ',a2)
31 format(4x,'| DFT exchange = ',f25.15,' ',a2)
32 format(4x,'| DFT correlation = ',f25.15,' ',a2)
33 format(4x,'| EXX contribution = ',f25.15,' ',a2)
40 format(4x,' |* pseudopotential energy = ',f25.15,' ',a2)
41 format(4x,' | core correction = ',f25.15,' ',a2)
42 format(4x,' | local contribution = ',f25.15,' ',a2)
43 format(4x,' | nonlocal contribution = ',f25.15,' ',a2)
60 format(4x,' |* pseudo/NA energy = ',f25.15,' ',a2)
62 format(4x,' | NA contribution = ',f25.15,' ',a2)
40 format(4x,'|* pseudopotential energy = ',f25.15,' ',a2)
41 format(4x,'| core correction = ',f25.15,' ',a2)
42 format(4x,'| local contribution = ',f25.15,' ',a2)
43 format(4x,'| nonlocal contribution = ',f25.15,' ',a2)
60 format(4x,'|* pseudo/NA energy = ',f25.15,' ',a2)
62 format(4x,'| NA contribution = ',f25.15,' ',a2)
11 format(4x,' | Ha correction = ',f25.15,' ',a2)
12 format(4x,' | XC correction = ',f25.15,' ',a2)
11 format(4x,'| Ha correction = ',f25.15,' ',a2)
12 format(4x,'| XC correction = ',f25.15,' ',a2)
10 format(4x,' |* Harris-Foulkes energy = ',f25.15,' ',a2)
13 format(4x,' |* DFT total energy = ',f25.15,' ',a2)
22 format(4x,' | estimated accuracy = ',f25.15,' ',a2)
10 format(4x,'|* Harris-Foulkes energy = ',f25.15,' ',a2)
13 format(4x,'|* DFT total energy = ',f25.15,' ',a2)
22 format(4x,'| estimated accuracy = ',f25.15,' ',a2)
14 format(4x,' | GS Energy as E-(1/2)TS = ',f25.15,' ',a2)
15 format(4x,' | Free Energy as E-TS = ',f25.15,' ',a2)
16 format(4x,' | GS Energy with kT -> 0 = ',f25.15,' ',a2)
17 format(4x,' | Dispersion (DFT-D2) = ',f25.15,' ',a2)
18 format(4x,' | cDFT Energy as 2Tr[K.W] = ',f25.15,' ',a2)
19 format(4x,' | Number of e- spin up = ',f25.15)
20 format(4x,' | Number of e- spin down = ',f25.15)
21 format(4x,' | Spin pol. as (up - down) = ',f25.15)
14 format(4x,'| GS Energy as E-(1/2)TS = ',f25.15,' ',a2)
15 format(4x,'| Free Energy as E-TS = ',f25.15,' ',a2)
16 format(4x,'| GS Energy with kT -> 0 = ',f25.15,' ',a2)
17 format(4x,'| Dispersion (DFT-D2) = ',f25.15,' ',a2)
18 format(4x,'| cDFT Energy as 2Tr[K.W] = ',f25.15,' ',a2)
19 format(4x,'| Number of e- spin up = ',f25.15)
20 format(4x,'| Number of e- spin down = ',f25.15)
21 format(4x,'| Spin pol. as (up - down) = ',f25.15)
23 format(4x,' |* check for accuracy = ',f25.15,' ',a2)
24 format(4x,' | number of e- num. int. = ',f25.15,' ',a2)
25 format(4x,' | number of e- as 2Tr[KS] = ',f25.15,' ',a2)
26 format(4x,' | one-electron energy = ',f25.15,' ',a2)
27 format(4x,' | potential energy V = ',f25.15,' ',a2)
28 format(4x,' | kinetic energy T = ',f25.15,' ',a2)
29 format(4x,' |* virial V/T = ',f25.15,' ',a2)
50 format(4x,' | rescaled DFT exchange = ',f25.15,' ',a2)
51 format(4x,' | rescaled exact exchange = ',f25.15,' ',a2)
23 format(4x,'|* check for accuracy = ',f25.15,' ',a2)
24 format(4x,'| number of e- num. int. = ',f25.15,' ',a2)
25 format(4x,'| number of e- as 2Tr[KS] = ',f25.15,' ',a2)
26 format(4x,'| one-electron energy = ',f25.15,' ',a2)
27 format(4x,'| potential energy V = ',f25.15,' ',a2)
28 format(4x,'| kinetic energy T = ',f25.15,' ',a2)
29 format(4x,'|* virial V/T = ',f25.15,' ',a2)
50 format(4x,'| rescaled DFT exchange = ',f25.15,' ',a2)
51 format(4x,'| rescaled exact exchange = ',f25.15,' ',a2)
end subroutine final_energy

148
src/force_module.f90
View File

@ -246,7 +246,7 @@ contains
flag_neutral_atom, flag_stress, &
rcellx, rcelly, rcellz, &
flag_atomic_stress, non_atomic_stress, &
flag_heat_flux, cell_constraint_flag
flag_heat_flux, cell_constraint_flag, min_layer
use density_module, only: get_electronic_density, density, &
build_Becke_weight_forces
use functions_on_grid, only: atomfns, H_on_atomfns
@ -261,8 +261,7 @@ contains
delta_E_xc, xc_energy, hartree_energy_drho
use hartree_module, only: Hartree_stress
use XC, ONLY: XC_GGA_stress
use input_module, only: leqi
use io_module, only: atom_output_threshold
use io_module, only: atom_output_threshold, return_prefix
implicit none
@ -292,7 +291,10 @@ contains
real(double), dimension(:), allocatable :: density_out_tot
real(double), dimension(:,:), allocatable :: density_out
character(len=1), dimension(3) :: comptstr = (/"x", "y", "z"/)
character(len=10) :: subname = "force: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
!****lat<$
if ( present(level) ) backtrace_level = level+1
if ( .not. present(level) ) backtrace_level = -10
@ -327,7 +329,7 @@ contains
call reg_alloc_mem(area_moveatoms, 3 * ni_in_cell, type_dbl)
cdft_force = zero
end if
if(iprint_MD>3) then
if(iprint_MD + min_layer>3) then
allocate(s_pulay_for(3,ni_in_cell),phi_pulay_for(3,ni_in_cell))
s_pulay_for = zero
phi_pulay_for = zero
@ -506,10 +508,11 @@ contains
max_force = zero
max_atom = 0
max_compt = 0
if (inode == ionode .and. write_forces .and. (iprint_MD>0 .or. ni_in_cell<atom_output_threshold)) then
write (io_lun, fmt='(/,4x,"Forces on atoms (",a2,"/",a2,")"/)') &
en_units(energy_units), d_units(dist_units)
write (io_lun, fmt='(4x," Atom X Y Z")')
if (inode == ionode .and. write_forces .and. (iprint_MD + min_layer>0 .or. ni_in_cell<atom_output_threshold)) then
write (io_lun, fmt='(/4x,a,a2,"/",a2,")"/)') &
trim(prefix)//" Forces on atoms (",en_units(energy_units), d_units(dist_units)
write (io_lun, fmt='(4x,a)') &
trim(prefix)//" Atom X Y Z"
end if
! Calculate forces and write out
g0 = zero
@ -550,44 +553,49 @@ contains
end if
end do ! j
if (inode == ionode) then
if(iprint_MD > 2) then
write(io_lun, 101) i
write(io_lun, 102) (for_conv * HF_force(j,i), j = 1, 3)
if(flag_neutral_atom_projector) write (io_lun, fmt='("Force NA : ",3f15.10)') (for_conv*NA_force(j,i),j=1,3)
write(io_lun, 112) (for_conv * HF_NL_force(j,i), j = 1, 3)
write(io_lun, 103) (for_conv * p_force(j,i), j = 1, 3)
if(iprint_MD>3) then
write (io_lun, fmt='(" Phi pulay : ",3f15.10)') (for_conv*phi_pulay_for(j,i),j=1,3)
write (io_lun, fmt='(" S pulay : ",3f15.10)') (for_conv*s_pulay_for(j,i),j=1,3)
if(iprint_MD + min_layer > 2) then
write(io_lun, 101) trim(prefix),i
write(io_lun, 102) trim(prefix),(for_conv * HF_force(j,i), j = 1, 3)
if(flag_neutral_atom_projector) &
write (io_lun, fmt='(4x,a,3f15.10)') trim(prefix)//"Force NA : ", &
(for_conv*NA_force(j,i),j=1,3)
write(io_lun, 112) trim(prefix),(for_conv * HF_NL_force(j,i), j = 1, 3)
write(io_lun, 103) trim(prefix),(for_conv * p_force(j,i), j = 1, 3)
if(iprint_MD + min_layer>3) then
write (io_lun, fmt='(4x,a,3f15.10)') trim(prefix)//" Phi pulay : ", &
(for_conv*phi_pulay_for(j,i),j=1,3)
write (io_lun, fmt='(4x,a,3f15.10)') trim(prefix)//" S pulay : ", &
(for_conv*s_pulay_for(j,i),j=1,3)
end if
write(io_lun, 104) (for_conv * KE_force(j,i), j = 1, 3)
write(io_lun, 104) trim(prefix),(for_conv * KE_force(j,i), j = 1, 3)
if(flag_neutral_atom) then
write(io_lun, 106) (for_conv * screened_ion_force(j,i), j = 1, 3)
write(io_lun, 106) trim(prefix),(for_conv * screened_ion_force(j,i), j = 1, 3)
else
write(io_lun, 106) (for_conv * ion_interaction_force(j,i), j = 1, 3)
write(io_lun, 106) trim(prefix),(for_conv * ion_interaction_force(j,i), j = 1, 3)
end if
if (flag_pcc_global) write(io_lun, 108) (for_conv * pcc_force(j,i), j = 1, 3)
if (flag_dft_d2) write (io_lun, 109) (for_conv * disp_force(j,i), j = 1, 3)
if (flag_perform_cdft) write (io_lun, fmt='("Force cDFT : ",3f15.10)') &
(for_conv*cdft_force(j,i),j=1,3)
if (flag_pcc_global) write(io_lun, 108) trim(prefix),(for_conv * pcc_force(j,i), j = 1, 3)
if (flag_dft_d2) write (io_lun, 109) trim(prefix),(for_conv * disp_force(j,i), j = 1, 3)
if (flag_perform_cdft) write (io_lun, fmt='(4x,a,3f15.10)') &
trim(prefix)//"Force cDFT : ",(for_conv*cdft_force(j,i),j=1,3)
if (flag_self_consistent) then
write (io_lun, 105) (for_conv * tot_force(j,i), j = 1, 3)
write (io_lun, 105) trim(prefix),(for_conv * tot_force(j,i), j = 1, 3)
else
write (io_lun, 107) (for_conv * nonSC_force(j,i), j = 1, 3)
write (io_lun, 105) (for_conv * tot_force(j,i), j = 1, 3)
write (io_lun, 107) trim(prefix),(for_conv * nonSC_force(j,i), j = 1, 3)
write (io_lun, 105) trim(prefix),(for_conv * tot_force(j,i), j = 1, 3)
end if
else if (write_forces .and. iprint_MD>0 .or. ni_in_cell<atom_output_threshold) then
write (io_lun,fmt='(6x,i6,3f15.10)') &
i, (for_conv * tot_force(j,i), j = 1, 3)
end if ! (iprint_MD > 2)
else if (write_forces .and. iprint_MD + min_layer>0 .or. ni_in_cell<atom_output_threshold) then
write (io_lun,fmt='(4x,a,i6,3f15.10)') &
trim(prefix),i, (for_conv * tot_force(j,i), j = 1, 3)
end if ! (iprint_MD + min_layer > 2)
end if ! (inode == ionode)
end do ! i
if (inode == ionode) then
write (io_lun, fmt='(4x,"Maximum force : ",f15.8,"(",a2,"/",&
write (io_lun, fmt='(4x,a,f15.8,"(",a2,"/",&
&a2,") on atom ",i9," in ",a1," direction")') &
for_conv * max_force, en_units(energy_units), &
trim(prefix)//" Maximum force : ",for_conv * max_force, en_units(energy_units), &
d_units(dist_units), max_atom, comptstr(max_compt)
write(io_lun, fmt='(4x,"Force Residual: ",f20.10," ",a2,"/",a2)') &
write(io_lun, fmt='(4x,a,f20.10," ",a2,"/",a2)') &
trim(prefix)//" Force Residual: ", &
for_conv*sqrt(g0/ni_in_cell), en_units(energy_units), d_units(dist_units)
end if
! We will add PCC and nonSCF stresses even if the flags are not set, as they are
@ -655,42 +663,42 @@ contains
end if
! Output
if (inode == ionode) then
write (io_lun,fmt='(/4x," ",3a15)') "X","Y","Z"
if(iprint_MD > 1) write(io_lun,fmt='(4x,"Stress contributions:")')
write (io_lun,fmt='(/4x,a,3a15)') trim(prefix)//" ","X","Y","Z"
if(iprint_MD + min_layer > 1) write(io_lun,fmt='(4x,a)') trim(prefix)//"Stress contributions:"
end if
call print_stress("K.E. stress: ", KE_stress, 3)
call print_stress("S-Pulay stress: ", SP_stress, 3)
call print_stress("Phi-Pulay stress: ", PP_stress, 3)
call print_stress("Local stress: ", loc_HF_stress, 3)
call print_stress("Local G stress: ", loc_G_stress, 3)
call print_stress("Non-local stress: ", NL_stress, 3)
call print_stress("Jacobian stress: ", GPV_stress, 3)
call print_stress("XC stress: ", XC_stress, 3)
call print_stress(trim(prefix)//" K.E. stress: ", KE_stress, 3 + min_layer)
call print_stress(trim(prefix)//" S-Pulay stress: ", SP_stress, 3 + min_layer)
call print_stress(trim(prefix)//" Phi-Pulay stress: ", PP_stress, 3 + min_layer)
call print_stress(trim(prefix)//" Local stress: ", loc_HF_stress, 3 + min_layer)
call print_stress(trim(prefix)//" Local G stress: ", loc_G_stress, 3 + min_layer)
call print_stress(trim(prefix)//" Non-local stress: ", NL_stress, 3 + min_layer)
call print_stress(trim(prefix)//" Jacobian stress: ", GPV_stress, 3 + min_layer)
call print_stress(trim(prefix)//" XC stress: ", XC_stress, 3 + min_layer)
if (flag_neutral_atom) then
call print_stress("Ion-ion stress: ", screened_ion_stress, 3)
call print_stress("N.A. stress: ", NA_stress, 3)
call print_stress(trim(prefix)//" Ion-ion stress: ", screened_ion_stress, 3 + min_layer)
if(flag_neutral_atom_projector) &
call print_stress(trim(prefix)//" N.A. stress: ", NA_stress, 3 + min_layer)
else
call print_stress("Ion-ion stress: ", ion_interaction_stress, 3)
call print_stress(trim(prefix)//" Ion-ion stress: ", ion_interaction_stress, 3 + min_layer)
end if
call print_stress("Hartree stress: ", Hartree_stress, 3)
call print_stress("PCC stress: ", pcc_stress, 3)
call print_stress("non-SCF stress: ", nonSCF_stress, 3)
if(flag_dft_D2) call print_stress("DFT-D2 stress: ", disp_stress, 3)
call print_stress("Total stress: ", stress, 0)
call print_stress(trim(prefix)//" Hartree stress: ", Hartree_stress, 3 + min_layer)
call print_stress(trim(prefix)//" PCC stress: ", pcc_stress, 3 + min_layer)
call print_stress(trim(prefix)//" non-SCF stress: ", nonSCF_stress, 3 + min_layer)
call print_stress(trim(prefix)//" Total stress: ", stress, 0 + min_layer)
volume = rcellx*rcelly*rcellz
! We need pressure in GPa, and only diagonal terms output
scale = -HaBohr3ToGPa/volume
!call print_stress("Total pressure: ", stress*scale, 0)
if(inode==ionode.AND.iprint_MD>=0) &
write(io_lun,'(/4x,a18,f15.8,a4)') "Average pressure: ", &
if(inode==ionode.AND.iprint_MD + min_layer>=0) &
write(io_lun,'(/4x,a,f15.8,a4)') trim(prefix)//" Average pressure: ", &
third*scale*(stress(1,1) + stress(2,2) + stress(3,3))," GPa"
if (flag_atomic_stress .and. iprint_MD > 2) call check_atomic_stress
if (flag_atomic_stress .and. iprint_MD + min_layer > 2) call check_atomic_stress
end if
call my_barrier()
if(iprint_MD>3) deallocate(s_pulay_for,phi_pulay_for)
if (inode == ionode .and. iprint_MD > 1 .and. write_forces) &
write (io_lun, fmt='(4x,"Finished force")')
if(iprint_MD + min_layer>3) deallocate(s_pulay_for,phi_pulay_for)
if (inode == ionode .and. iprint_MD + min_layer > 1 .and. write_forces) &
write (io_lun, fmt='(4x,a)') trim(prefix)//"Finished force"
call start_timer(tmr_std_allocation)
if (flag_pcc_global) then
@ -717,16 +725,16 @@ contains
return
101 format('Force on atom ',i9)
102 format('Force H-F : ',3f15.10)
112 format('Force H-Fnl : ',3f15.10)
103 format('Force pulay : ',3f15.10)
104 format('Force KE : ',3f15.10)
106 format('Force Ion-Ion: ',3f15.10)
107 format('Force nonSC : ',3f15.10)
108 format('Force PCC : ',3f15.10)
105 format('Force Total : ',3f15.10)
109 format('Force disp : ',3f15.10)
101 format(a,'Force on atom ',i9)
102 format(a,'Force H-F : ',3f15.10)
112 format(a,'Force H-Fnl : ',3f15.10)
103 format(a,'Force pulay : ',3f15.10)
104 format(a,'Force KE : ',3f15.10)
106 format(a,'Force Ion-Ion: ',3f15.10)
107 format(a,'Force nonSC : ',3f15.10)
108 format(a,'Force PCC : ',3f15.10)
105 format(a,'Force Total : ',3f15.10)
109 format(a,'Force disp : ',3f15.10)
end subroutine force
!!***
@ -4312,7 +4320,7 @@ subroutine print_stress(label, str_mat, print_level)
integer, intent(in) :: print_level
! local variables
character(20) :: fmt = '(4x,a18,3f15.8,a3)'
character(20) :: fmt = '(4x,a,3f15.8,a3)'
character(20) :: blank = ''
if (inode==ionode) then

1
src/global_module.f90
View File

@ -293,6 +293,7 @@ module global_module
integer, parameter :: IPRINT_TIME_THRES2 = 4 ! Not that important
integer, parameter :: IPRINT_TIME_THRES3 = 6 ! For special purposes
integer :: min_layer ! Layer of minimisation algorithm (from 0 to -n)
! For P.C.C.
logical :: flag_pcc_global = .false.

3
src/initial_read_module.f90
View File

@ -1325,6 +1325,7 @@ contains
charge_up(i) = fdf_double ('Atom.SpinNeUp',zero)
charge_dn(i) = fdf_double ('Atom.SpinNeDn',zero)
sum_elecN_spin = charge_up(i)+charge_dn(i)
flag_InitialAtomicSpin = .false.
if (abs(sum_elecN_spin)>RD_ERR) then
flag_InitialAtomicSpin = .true.
! We will check that the sum of charge_up and charge_dn matches charge later
@ -2635,7 +2636,7 @@ contains
write(io_lun,fmt='(10x,"SCF tolerance: ",f12.8)') sc_tolerance
write(io_lun,fmt='(10x,"SCF iterations: ",i4)') maxitersSC
end if
if(flag_self_consistent) then
if(.not.flag_diagonalisation) then
write(io_lun,fmt='(10x,"O(N) tolerance: ",f12.8)') L_tolerance
write(io_lun,fmt='(10x,"O(N) iterations: ",i4)') n_L_iterations
end if

120
src/initialisation_module.f90
View File

@ -122,7 +122,7 @@ contains
flag_only_dispersion, flag_neutral_atom, &
flag_atomic_stress, flag_heat_flux, &
flag_full_stress, area_moveatoms, &
atomic_stress, non_atomic_stress
atomic_stress, non_atomic_stress, min_layer
use GenComms, only: inode, ionode, my_barrier, end_comms, &
cq_abort
use initial_read, only: read_and_write
@ -137,6 +137,7 @@ contains
use angular_coeff_routines, only: set_fact
use maxima_module, only: lmax_ps, lmax_pao
use XC, only: init_xc
use io_module, only: return_prefix
implicit none
@ -151,7 +152,10 @@ contains
logical :: start, start_L
logical :: read_phi
integer :: lmax_tot, stat
character(len=12) :: subname = "initialise: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
call init_timing_system(inode)
call init_reg_mem
@ -199,13 +203,13 @@ contains
if (flag_heat_flux) then
if (.not. flag_full_stress) then
flag_full_stress = .true.
if (inode==ionode) write(io_lun,'(2x,a)') &
"WARNING: setting AtomMove.FullStress T for heat flux calculation"
if (inode==ionode) write(io_lun,'(4x,a)') &
trim(prefix)//"WARNING: setting AtomMove.FullStress T for heat flux calculation"
end if
if (.not. flag_atomic_stress) then
flag_atomic_stress = .true.
if (inode==ionode) write(io_lun,'(2x,a)') &
"WARNING: setting AtomMove.AtomicStress T for heat flux calculation"
if (inode==ionode) write(io_lun,'(4x,a)') &
trim(prefix)//"WARNING: setting AtomMove.AtomicStress T for heat flux calculation"
end if
end if
if (flag_atomic_stress) then
@ -323,7 +327,7 @@ contains
flag_Becke_weights, &
flag_pcc_global, flag_dft_d2, &
iprint_gen, flag_perform_cDFT, &
nspin, &
nspin, min_layer, &
glob2node, flag_XLBOMD, &
flag_neutral_atom, flag_diagonalisation
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, &
@ -364,7 +368,7 @@ contains
grid_point_position
use primary_module, only: bundle
use group_module, only: blocks
use io_module, only: read_blocks
use io_module, only: read_blocks, return_prefix
use functions_on_grid, only: associate_fn_on_grid
use potential_module, only: potential
use maxima_module, only: maxngrid, lmax_ps, lmax_pao
@ -390,6 +394,8 @@ contains
real(double) :: rcut_BCS !TM 26/Jun/2003
type(cq_timer) :: backtrace_timer
integer :: backtrace_level
character(len=7) :: subname = "setup: "
character(len=120) :: prefix
!****lat<$
if ( present(level) ) backtrace_level = level+1
@ -397,6 +403,7 @@ contains
call start_backtrace(t=backtrace_timer,who='set_up',&
where=area,level=backtrace_level,echo=.true.)
!****lat>$
prefix = return_prefix(subname, min_layer)
! Set organisation of blocks of grid-points.
! set_blocks determines the number of blocks on this node,
@ -409,8 +416,8 @@ contains
call set_blocks_from_new()
! Allocate ?
!call set_blocks(inode, ionode)
if (inode == ionode .and. iprint_init > 1) &
write(io_lun,*) 'Completed set_blocks()'
if (inode == ionode .and. iprint_init + min_layer > 2) &
write(io_lun,fmt='(4x,a)') trim(prefix)//'Completed set_blocks()'
n_my_grid_points = n_blocks*n_pts_in_block
! allocate(grid_point_x(n_my_grid_points),&
! grid_point_y(n_my_grid_points),&
@ -450,14 +457,14 @@ contains
call reg_alloc_mem(area_index, maxngrid, type_dbl)
! extra local potential for second spin channel for spin polarised calculation
call my_barrier()
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed set_domains()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed set_domains()'
! Sorts out which processor owns which atoms
call distribute_atoms(inode, ionode)
call my_barrier
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed distribute_atoms()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed distribute_atoms()'
! Create a covering set
call my_barrier
!Define rcut_BCS !TM 26/Jun/2003
@ -466,8 +473,8 @@ contains
! if(rcut_BCS < rcut(i)) rcut_BCS= rcut(i)
!enddo ! i=1, mx_matrices
rcut_BCS = rcut(max_range)
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) ' rcut for BCS_parts =', rcut_BCS
if (inode == ionode .and. iprint_init > 3) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'rcut for BCS_parts =', rcut_BCS
call make_cs(inode-1, rcut_BCS, BCS_parts, parts, bundle, &
ni_in_cell, x_atom_cell, y_atom_cell, z_atom_cell)
@ -475,16 +482,16 @@ contains
call make_iprim(BCS_parts, bundle)
call send_ncover(BCS_parts, inode)
call my_barrier
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Made covering set for matrix multiplications'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Made covering set for matrix multiplications'
! Create all of the indexing required to perform matrix multiplications
! at a later point. This routine also identifies all the density
! matrix range interactions and hamiltonian range interactions
! associated with any atom being handled by this processor.
call immi(parts, bundle, BCS_parts, inode)
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed immi()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed immi()'
if (flag_XLBOMD) call immi_XL(parts,bundle,BCS_parts,inode)
! set up all the data block by block for atoms overlapping any
@ -495,13 +502,13 @@ contains
call setgrid(inode-1, r_core_squared, r_h)
call my_barrier()
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed set_grid()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed set_grid()'
call associate_fn_on_grid
call my_barrier()
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed associate_fn_on_grid()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed associate_fn_on_grid()'
! The FFT requires the data to be reorganised into columns parallel to
! each axis in turn. The data for this organisation is help in map.inc,
@ -524,34 +531,34 @@ contains
call cq_abort("Error deallocating chdenr: ", maxngrid, stat)
call reg_dealloc_mem(area_init, maxngrid, type_dbl)
call my_barrier()
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed fft init'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed fft init'
! Initialise the routines to calculate ion-ion interactions
if(flag_neutral_atom) then
call setup_screened_ion_interaction
call my_barrier
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed setup_ion_interaction()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed setup_ion_interaction()'
else
! set up the Ewald sumation: find out how many superlatices
! in the real space sum and how many reciprocal latice vectors in the
! reciprocal space sum are needed for a given energy tolerance.
call set_ewald(inode,ionode)
call my_barrier
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Completed set_ewald()'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Completed set_ewald()'
end if
! +++
! Generate D2CS
if (flag_dft_d2) then
if ((inode == ionode) .and. (iprint_gen > 1) ) &
if ((inode == ionode) .and. (iprint_gen > 2) ) &
write (io_lun, '(/1x,"The dispersion is considered in the &
&DFT-D2 level.")')
call make_cs(inode-1, r_dft_d2, D2_CS, parts, bundle, ni_in_cell, &
x_atom_cell, y_atom_cell, z_atom_cell)
if ( (inode == ionode) .and. (iprint_gen > 1) ) then
if ( (inode == ionode) .and. (iprint_gen > 3) ) then
write (io_lun, '(/8x,"+++ D2_CS%ng_cover:",i10)') &
D2_CS%ng_cover
write (io_lun, '(8x,"+++ D2_CS%ncoverx, y, z:",3i8)') &
@ -562,11 +569,6 @@ contains
D2_CS%nx_origin, D2_CS%ny_origin, D2_CS%nz_origin
end if
call set_para_D2
if (inode == ionode) then !! DEBUG !!
write (io_lun, '(a, f10.5)') & !! DEBUG !!
"Sbrt: make_cs for DFT-D2, the cutoff is ", & !! DEBUG !!
r_dft_d2 !! DEBUG !!
end if !! DEBUG !!
end if
! external potential - first set up angular momentum bits
@ -621,8 +623,8 @@ contains
call build_Becke_weights
call build_Becke_charges(atomcharge, density, maxngrid)
end if
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Done init_pseudo '
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') trim(prefix)//'Done init_pseudo '
if(flag_diagonalisation) then
call init_blacs_pg
@ -1226,7 +1228,7 @@ contains
! make SF-PAO coefficients
call initial_SFcoeff(.true., .true., fixed_potential, .true.)
endif
if (inode == ionode .and. iprint_init > 1) write (io_lun, *) 'Got SFcoeff'
if (inode == ionode .and. iprint_init > 2) write (io_lun, fmt='(4x,a)') 'Got SFcoeff'
call my_barrier
endif
!!$
@ -1248,7 +1250,7 @@ contains
else
call get_S_matrix(inode, ionode)
endif
if (inode == ionode .and. iprint_init > 1) write (io_lun, *) 'Got S'
if (inode == ionode .and. iprint_init > 2) write (io_lun, fmt='(4x,a)') 'Got S'
call my_barrier
!!$
!!$
@ -1260,8 +1262,8 @@ contains
if (.not. flag_diagonalisation .and. find_chdens .and. (start .or. start_L)) then
call initial_L()
call my_barrier()
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Got L matrix'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun, fmt='(4x,a)') 'Got L matrix'
if (vary_mu) then
! This cannot be timed within the routine
call start_timer(tmr_std_densitymat)
@ -1274,12 +1276,12 @@ contains
call grab_matrix2('L',inode,nfile,Info,InfoGlob,index=index_MatrixFile,n_matrix=nspin)
call my_barrier()
call Matrix_CommRebuild(InfoGlob,Info,Lrange,L_trans,matL,nfile,symm,n_matrix=nspin)
if (inode == ionode .and. iprint_init > 1) write (io_lun, *) 'Grabbed L matrix'
if (inode == ionode .and. iprint_init > 2) write (io_lun, fmt='(4x,a)') 'Grabbed L matrix'
else
call grab_matrix2('K',inode,nfile,Info,InfoGlob,index=index_MatrixFile,n_matrix=nspin)
call my_barrier()
call Matrix_CommRebuild(InfoGlob,Info,Hrange,H_trans,matK,nfile,n_matrix=nspin)
if (inode == ionode .and. iprint_init > 1) write (io_lun, *) 'Grabbed K matrix'
if (inode == ionode .and. iprint_init > 2) write (io_lun, fmt='(4x,a)') 'Grabbed K matrix'
!DEBUG call Report_UpdateMatrix("Kmat")
end if
end if
@ -1303,8 +1305,8 @@ contains
dontM2, dontM3, dontM4, dophi, dontE, &
mat_phi=matphi)
electrons_tot = spin_factor * sum(electrons)
if (inode == ionode .and. iprint_init > 1) &
write (io_lun,*) 'Got elect: (Nup, Ndn, Ntotal) ', &
if (inode == ionode .and. iprint_init > 2) &
write (io_lun,fmt='(4x,a)') 'Got elect: (Nup, Ndn, Ntotal) ', &
electrons(1), electrons(nspin), &
electrons_tot
end if
@ -1319,15 +1321,15 @@ contains
!!$
!!$
! (5) Find the Ewald energy for the initial set of atoms
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Ionic electrostatics'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun, fmt='(4x,a)') 'Ionic electrostatics'
if(flag_neutral_atom) then
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Calling screened_ion_interaction'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun, fmt='(4x,a)') 'Calling screened_ion_interaction'
call screened_ion_interaction
else
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'Calling ewald'
if (inode == ionode .and. iprint_init > 2) &
write (io_lun, fmt='(4x,a)') 'Calling ewald'
call ewald
end if
!!$
@ -1337,8 +1339,8 @@ contains
! +++
! (6) Find the dispersion energy for the initial set of atoms
if (flag_dft_d2) then
if ((inode == ionode) .and. (iprint_init > 1) ) &
write (io_lun, *) 'Calling DFT-D2'
if ((inode == ionode) .and. (iprint_init > 2) ) &
write (io_lun, fmt='(4x,a)') 'Calling DFT-D2'
call dispersion_D2
end if
call my_barrier
@ -1347,7 +1349,7 @@ contains
!!$
!!$
if (inode == ionode .and. iprint_init > 2) &
write (io_lun, *) 'Find_chdens is ', find_chdens
write (io_lun, fmt='(4x,a)') 'Find_chdens is ', find_chdens
!!$
!!$
!!$
@ -1359,8 +1361,8 @@ contains
maxngrid)
electrons_tot = spin_factor * sum(electrons)
density = density * ne_in_cell/electrons_tot
if (inode == ionode .and. iprint_init > 1) &
write (io_lun, *) 'In initial_H, electrons: ', electrons_tot
if (inode == ionode .and. iprint_init > 2) &
write (io_lun, fmt='(4x,a)') 'In initial_H, electrons: ', electrons_tot
! if flag_LFD=T, update SF-PAO coefficients with the obtained density unless they have been read
! and update S with the coefficients
if ((.NOT.read_option).AND.flag_LFD) then
@ -1557,8 +1559,8 @@ contains
real(double) :: rcut_max, r_core
!-- Start of the subroutine (set_grid_new)
if(myid == 0 .and. iprint_index > 1) &
write (io_lun, *) 'setgrid_new starts'
if(myid == 0 .and. iprint_index > 2) &
write (io_lun, fmt='(4x,a)') 'setgrid_new starts'
!if(iprint_index > 4) write(io_lun,*) ' setgrid_new starts for myid= ',myid
!Sets up domain

12
src/io_module.f90
View File

@ -812,7 +812,7 @@ second: do
!call create_sfc_partitions(myid, parts)
call sfc_partitions_to_processors(parts)
np_in_cell = parts%ngcellx*parts%ngcelly*parts%ngcellz
if (iprint_init > 1.AND.myid==0) write(io_lun,*) 'Finished partitioning'
if (iprint_init > 2.AND.myid==0) write(io_lun,fmt='(4x,a)') 'Finished partitioning'
end if
! inverse table to npnode
do np=1,np_in_cell
@ -3453,5 +3453,15 @@ second: do
end subroutine check_stop
!!***
function return_prefix(name, level)
character(len=120) :: return_prefix
character(len=*) :: name
character(len=10):: prefix = " "
integer :: level
return_prefix = prefix(1:-2*level)//name
end function return_prefix
end module io_module

10
src/minimise.f90
View File

@ -149,7 +149,7 @@ contains
flag_LmatrixReuse,McWFreq, &
flag_multisite, &
io_lun, flag_out_wf, wf_self_con, flag_write_DOS, &
flag_diagonalisation, nspin, flag_LFD
flag_diagonalisation, nspin, flag_LFD, min_layer
use energy, only: get_energy, xc_energy, final_energy
use GenComms, only: cq_abort, inode, ionode
use blip_minimisation, only: vary_support, dE_blip
@ -237,7 +237,7 @@ contains
call vary_pao(n_support_iterations, fixed_potential, &
vary_mu, n_L_iterations, L_tolerance, &
sc_tolerance, energy_tolerance, &
total_energy, expected_reduction, level)
total_energy, expected_reduction)
end if
dE_elec_opt = dE_PAO
else ! Or SCF if necessary
@ -255,7 +255,7 @@ contains
dE_elec_opt = dE_blip
else if (flag_basis_set == PAOs) then
if (flag_multisite .and. flag_LFD_NonSCF) then
if (inode==ionode) write(io_lun,'(/3x,A/)') &
if (inode==ionode) write(io_lun,'(/4x,A/)') &
'WARNING: Numerical PAO minimisation will be performed without doing LFD_SCF!'
!2017.Dec.28 TM: We need Selfconsistent Hamiltonian if this routine is called from control
call new_SC_potl(.false., sc_tolerance, reset_L, &
@ -272,7 +272,7 @@ contains
call vary_pao(n_support_iterations, fixed_potential, &
vary_mu, n_L_iterations, L_tolerance, &
sc_tolerance, energy_tolerance, &
total_energy, expected_reduction, level)
total_energy, expected_reduction)
end if
dE_elec_opt = dE_PAO
else
@ -305,7 +305,7 @@ contains
call vary_pao(n_support_iterations, fixed_potential, &
vary_mu, n_L_iterations, L_tolerance, &
sc_tolerance, energy_tolerance, &
total_energy, expected_reduction, level)
total_energy, expected_reduction)
end if
else
call cq_abort("get_E_and_F: basis set undefined: ", &

105
src/move_atoms.module.f90
View File

@ -760,7 +760,7 @@ contains
IPRINT_TIME_THRES1, flag_pcc_global, &
id_glob, &
flag_LmatrixReuse, flag_diagonalisation, nspin, &
flag_SFcoeffReuse
flag_SFcoeffReuse, min_layer
use minimise, only: get_E_and_F, sc_tolerance, L_tolerance, &
n_L_iterations
use GenComms, only: my_barrier, myid, inode, ionode, &
@ -806,6 +806,7 @@ contains
integer :: ig, both, mat, update_var
write(io_lun,*) 'min_layer is ',min_layer
call start_timer(tmr_std_moveatoms)
if(flag_SFcoeffReuse) then
@ -876,8 +877,10 @@ contains
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
end if
min_layer = min_layer - 1
call get_E_and_F(fixed_potential, vary_mu, e3, .false., &
.false.)
min_layer = min_layer + 1
! Now, we call dump_pos_and_matrices here. : 2018.Jan19 TM
! but if we want to use the information of the matrices in the beginning of this line minimisation
! you can comment the following line, in the future.
@ -969,7 +972,8 @@ contains
IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground state before
! we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
@ -980,6 +984,7 @@ contains
else
call get_E_and_F(fixed_potential, vary_mu, energy_out, .true., .false.)
end if
min_layer = min_layer + 1
! 2018.Jan19 TM
call dump_pos_and_matrices
@ -1062,7 +1067,8 @@ contains
IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if(flag_vary_basis) then
min_layer = min_layer - 1
if(flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
@ -1075,7 +1081,7 @@ contains
call get_E_and_F(fixed_potential, vary_mu, energy_out, &
.true., .false.)
end if
min_layer = min_layer + 1
! 2018.Jan19 TM : probably we don't need to call dump_pos_and_matrices here, since
! we will call it after calling safemin2
call dump_pos_and_matrices
@ -1128,7 +1134,7 @@ contains
IPRINT_TIME_THRES1, flag_pcc_global, &
id_glob, &
flag_LmatrixReuse, flag_diagonalisation, nspin, &
flag_SFcoeffReuse
flag_SFcoeffReuse, min_layer
use minimise, only: get_E_and_F, sc_tolerance, L_tolerance, &
n_L_iterations
use GenComms, only: my_barrier, myid, inode, ionode, &
@ -1146,7 +1152,7 @@ contains
use store_matrix, ONLY: dump_pos_and_matrices
use mult_module, ONLY: allocate_temp_matrix, free_temp_matrix, matrix_sum
use global_module, ONLY: atomf, sf
use io_module, ONLY: dump_matrix
use io_module, ONLY: dump_matrix, return_prefix
use force_module, only: force
implicit none
@ -1173,15 +1179,19 @@ contains
integer :: ig, both, mat
character(len=10) :: subname = "back_lm: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
call start_timer(tmr_std_moveatoms)
iter = 0
old_alpha = zero
alpha = one
e0 = total_energy
if (inode == ionode .and. iprint_MD > 0) &
if (inode == ionode .and. iprint_MD + min_layer > 0) &
write (io_lun, &
fmt='(4x,"In backtrack_linemin, initial energy is ",f16.6," ",a2)') &
fmt='(4x,a,f16.6," ",a2)') trim(prefix)//" Initial energy is ",&
en_conv * energy_in, en_units(energy_units)
c1 = 0.1_double
@ -1194,9 +1204,9 @@ contains
grad_f_dot_p = grad_f_dot_p - direction(2,i)*tot_force(2,j)
grad_f_dot_p = grad_f_dot_p - direction(3,i)*tot_force(3,j)
end do
if(inode==ionode.AND.iprint_MD>1) &
write(io_lun, fmt='(2x,"Starting backtrack_linemin, magnitude of grad_f.p is ",e16.6)') &
sqrt(-grad_f_dot_p/ni_in_cell)
if(inode==ionode.AND.iprint_MD + min_layer>1) &
write(io_lun, fmt='(4x,a,e16.6)') &
trim(prefix)//" Magnitude of grad_f.p is ",sqrt(-grad_f_dot_p/ni_in_cell)
done = .false.
do while (.not. done)
iter = iter+1
@ -1213,37 +1223,40 @@ contains
else
call update_pos_and_matrices(updateLorK,direction)
endif
if (inode == ionode .and. iprint_MD > 2) then
if (inode == ionode .and. iprint_MD + min_layer > 2) then
do i=1,ni_in_cell
write (io_lun,fmt='(2x,"Position: ",i3,3f13.8)') i, &
write (io_lun,fmt='(4x,a,i3,3f13.8)') trim(prefix)//"Position: ",i, &
x_atom_cell(i), y_atom_cell(i), z_atom_cell(i)
end do
end if
call update_H(fixed_potential)
! Write out atomic positions
if (iprint_MD > 2) then
if (iprint_MD + min_layer > 2) then
call write_atomic_positions("UpdatedAtoms_tmp.dat", &
trim(pdb_template))
end if
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
end if
call get_E_and_F(fixed_potential, vary_mu, e3, .false., &
.false.)
min_layer = min_layer + 1
!call dump_pos_and_matrices
! e3 is f(x + alpha p)
armijo = e0 + c1 * alpha * grad_f_dot_p
if (inode == ionode .and. iprint_MD > 1) then
if (inode == ionode .and. iprint_MD + min_layer > 1) then
write (io_lun, &
fmt='(4x,"In backtrack_linemin, iter ",i3," step and energy &
&are ",2f16.6," ",a2)') &
fmt='(4x,a,i3," step and energy &
&are ",2f16.6," ",a2)') trim(prefix)//" Iter ",&
iter, alpha, en_conv * e3, en_units(energy_units)
write(io_lun, fmt='(6x,"Armijo threshold is ",f16.6," ",a2)') armijo, en_units(energy_units)
write(io_lun, fmt='(4x,a,f16.6," ",a2)') trim(prefix)//"Armijo threshold is ", &
armijo, en_units(energy_units)
end if
if(e3<armijo) then ! success
done = .true.
@ -1256,8 +1269,10 @@ contains
energy_out = e3
call dump_pos_and_matrices
! Now find forces
min_layer = min_layer - 1
call force(fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, sc_tolerance, energy_out, .true.)
min_layer = min_layer + 1
! Evaluate new grad f dot p
grad_fp_dot_p = zero
do i=1, ni_in_cell
@ -1266,19 +1281,17 @@ contains
grad_fp_dot_p = grad_f_dot_p - direction(2,i)*tot_force(2,j)
grad_fp_dot_p = grad_f_dot_p - direction(3,i)*tot_force(3,j)
end do
if(inode==ionode.AND.iprint_MD>3) &
write(io_lun,fmt='(6x,"In backtrack_linemin, second Wolfe: ",e11.4," < ",e11.4)') &
if(inode==ionode.AND.iprint_MD + min_layer>3) &
write(io_lun,fmt='(4x,a,e11.4," < ",e11.4)') &
trim(prefix)//" Second Wolfe: ",&
abs(grad_fp_dot_p), c2*abs(grad_f_dot_p)
dE = e0 - energy_out
if (inode == ionode .and. iprint_MD > 2) then
if (inode == ionode .and. iprint_MD + min_layer >= 0) then
write (io_lun, &
fmt='(4x,"In backtrack_linemin, exit after ",i4," &
&iterations with energy ",f16.6," ",a2)') &
fmt='(4x,a,i4," &
&iterations with energy ",f16.6," ",a2)') trim(prefix)//" Exit after ",&
iter, en_conv * energy_out, en_units(energy_units)
else if (inode == ionode .and. iprint_MD > 0) then
write (io_lun, fmt='(/4x,"In backtrack_linemin, final energy is ",f16.6," ",a2)') &
en_conv * energy_out, en_units(energy_units)
end if
call stop_timer(tmr_std_moveatoms)
@ -1470,7 +1483,7 @@ contains
IPRINT_TIME_THRES1, flag_pcc_global, &
id_glob, &
flag_LmatrixReuse, flag_diagonalisation, nspin, &
flag_SFcoeffReuse
flag_SFcoeffReuse, min_layer
use minimise, only: get_E_and_F, sc_tolerance, L_tolerance, &
n_L_iterations
use GenComms, only: my_barrier, myid, inode, ionode, &
@ -1569,13 +1582,15 @@ contains
end if
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
end if
call get_E_and_F(fixed_potential, vary_mu, e3, .false., &
.false.)
min_layer = min_layer + 1
!call dump_pos_and_matrices
! e3 is f(x + alpha p)
armijo = e0 + c1 * alpha * grad_f_dot_p
@ -1603,8 +1618,10 @@ contains
alpha = (one+scale)*alpha
end if
! Now find forces
min_layer = min_layer - 1
call force(fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, sc_tolerance, energy_out, .true.)
min_layer = min_layer + 1
grad_fp_dot_p = zero
do i=1, ni_in_cell
j = id_glob(i)
@ -1680,7 +1697,7 @@ contains
flag_reset_dens_on_atom_move, &
IPRINT_TIME_THRES1, flag_pcc_global, &
flag_diagonalisation, cell_constraint_flag, &
flag_SFcoeffReuse
flag_SFcoeffReuse, min_layer
use minimise, only: get_E_and_F, sc_tolerance, L_tolerance, &
n_L_iterations
use GenComms, only: my_barrier, myid, inode, ionode, cq_abort, &
@ -1781,13 +1798,15 @@ contains
call stop_print_timer(tmr_l_tmp1, "atom updates", IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
end if
call get_E_and_F(fixed_potential, vary_mu, e3, .false., &
.false.)
min_layer = min_layer + 1
call dump_pos_and_matrices
h3 = enthalpy(e3, target_press)
@ -1862,7 +1881,8 @@ contains
"safemin_cell - Final interpolation and updates", &
IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground state before we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
@ -1874,6 +1894,7 @@ contains
else
call get_E_and_F(fixed_potential, vary_mu, energy_out, .true., .false.)
end if
min_layer = min_layer + 1
call dump_pos_and_matrices
enthalpy_out = enthalpy(energy_out, target_press)
@ -1936,7 +1957,8 @@ contains
IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if(flag_vary_basis) then
min_layer = min_layer - 1
if(flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
@ -1948,6 +1970,7 @@ contains
call get_E_and_F(fixed_potential, vary_mu, energy_out, &
.true., .false.)
end if
min_layer = min_layer + 1
! we may not need to call dump_pos_and_matrices here. (if it would be called in the part after calling safemin_cell)
call dump_pos_and_matrices
enthalpy_out = enthalpy(energy_out, target_press)
@ -2004,7 +2027,7 @@ contains
flag_self_consistent, &
IPRINT_TIME_THRES1, flag_pcc_global, &
flag_LmatrixReuse, flag_SFcoeffReuse, &
atom_coord, id_glob
atom_coord, id_glob, min_layer
use minimise, only: get_E_and_F, sc_tolerance, L_tolerance, &
n_L_iterations
use GenComms, only: inode, ionode, cq_abort
@ -2122,13 +2145,15 @@ contains
call stop_print_timer(tmr_l_tmp1, "atom updates", IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
end if
call get_E_and_F(fixed_potential, vary_mu, e3, .false., &
.false.)
min_layer = min_layer + 1
! Now, we call dump_pos_and_matrices here. : 2018.Jan19 TM
! but if we want to use the information of the matrices in the beginning of this line minimisation
! you can comment the following line, in the future.
@ -2221,7 +2246,8 @@ contains
IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground state before
! we can minimise blips !
if (flag_vary_basis) then
min_layer = min_layer - 1
if (flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
@ -2232,6 +2258,7 @@ contains
else
call get_E_and_F(fixed_potential, vary_mu, energy_out, .true., .false.)
end if
min_layer = min_layer + 1
! 2018.Jan19 TM
call dump_pos_and_matrices
@ -2314,7 +2341,8 @@ contains
IPRINT_TIME_THRES1)
! We've just moved the atoms - we need a self-consistent ground
! state before we can minimise blips !
if(flag_vary_basis) then
min_layer = min_layer - 1
if(flag_vary_basis .or. .NOT.flag_LFD_NonSCF) then
call new_SC_potl(.false., sc_tolerance, reset_L, &
fixed_potential, vary_mu, n_L_iterations, &
L_tolerance, e3)
@ -2327,6 +2355,7 @@ contains
call get_E_and_F(fixed_potential, vary_mu, energy_out, &
.true., .false.)
end if
min_layer = min_layer + 1
! 2018.Jan19 TM : probably we don't need to call dump_pos_and_matrices here, since
! we will call it after calling safemin2
@ -4113,7 +4142,7 @@ contains
rcellz, flag_self_consistent, &
flag_reset_dens_on_atom_move, &
IPRINT_TIME_THRES1, flag_pcc_global, &
flag_diagonalisation, cell_constraint_flag
flag_diagonalisation, cell_constraint_flag, min_layer
use minimise, only: get_E_and_F, sc_tolerance, L_tolerance, &
n_L_iterations
use GenComms, only: my_barrier, myid, inode, ionode, &

149
src/multisiteSF_module.f90
View File

@ -45,7 +45,7 @@ module multisiteSF_module
use datatypes
use numbers
use global_module, only: io_lun, iprint_basis, area_basis
use global_module, only: io_lun, iprint_basis, area_basis, min_layer
use GenComms, only: myid, my_barrier, cq_abort, inode, ionode
use timer_module, only: start_timer,stop_timer
use timer_stdclocks_module, only: tmr_std_basis,tmr_std_matrices
@ -134,8 +134,8 @@ contains
use S_matrix_module, only: get_S_matrix
use H_matrix_module, only: get_H_matrix
use store_matrix, only: dump_pos_and_matrices
! use io_module, only: dump_matrix
use GenComms, only: mtime
use io_module, only: return_prefix
implicit none
@ -145,9 +145,12 @@ contains
real(double), dimension(nspin) :: electrons
real(double) :: electrons_tot, t0, t1, t2
integer :: spin_SF
character(len=14) :: subname = "Init_SFcoeff: "
character(len=120) :: prefix
if (inode==ionode .and. iprint_basis>3) write(io_lun,*) 'We are in sub:initial_SFcoeff'
prefix = return_prefix(subname, min_layer)
if (inode==ionode .and. iprint_basis + min_layer>2) write(io_lun,fmt='(4x,a)') &
trim(prefix)//' Entering'
if (output_naba_in_MSSF) call print_naba_in_MSSF
@ -169,8 +172,8 @@ contains
density, maxngrid, transform_AtomF_to_SF=.false.)
call my_barrier()
t1 = mtime()
if (inode == ionode .and. iprint_basis>4) write (io_lun,'(A,f20.8,A)') &
'LFD: Time for S and H in primitive PAO: ', t1-t0, ' ms'
if (inode == ionode .and. iprint_basis + min_layer>2) write (io_lun,'(4x,A,f20.8,A)') &
trim(prefix)//' Time for S and H in primitive PAO: ', t1-t0, ' ms'
t0 = t1
! Rayson's Localised Filter Diagonalisation method
call LocFilterDiag(mult(aLa_aHa_aLHa)%ahalo)
@ -200,12 +203,12 @@ contains
enddo
! Write out current SF coefficients with some iprint (in future)
if (iprint_basis>=3) call dump_pos_and_matrices(index=98)
if (iprint_basis + min_layer>=3) call dump_pos_and_matrices(index=98)
call my_barrier()
t2 = mtime()
if (inode == ionode .and. iprint_basis>2) write (io_lun,'(A,f20.8,A)') &
'Time for making initial multi-site SF coeffficients: ', t2-t0, ' ms'
if (inode == ionode .and. iprint_basis + min_layer>1) write (io_lun,'(4x,A,f20.8,A)') &
trim(prefix)//' Time for making initial multi-site SF coeffficients: ', t2-t0, ' ms'
return
end subroutine initial_SFcoeff
@ -258,7 +261,7 @@ contains
call start_timer(tmr_std_basis)
call start_timer(tmr_l_tmp1,WITH_LEVEL)
if(iprint_basis>=5.AND.myid==0) write(io_lun,'(6x,A)') ' We are in normalise_SFcoeff'
if(iprint_basis + min_layer>=5.AND.myid==0) write(io_lun,'(6x,A)') ' We are in normalise_SFcoeff'
call start_timer(tmr_std_matrices)
@ -396,27 +399,27 @@ contains
call start_timer(tmr_std_basis)
call start_timer(tmr_l_tmp1,WITH_LEVEL)
if(iprint_basis>=5.AND.myid==0) write(io_lun,'(6x,A)') ' We are in initial_SFcoeff_onsite'
if(iprint_basis + min_layer>=5.AND.myid==0) write(io_lun,'(6x,A)') ' We are in initial_SFcoeff_onsite'
if(iprint_basis>=5.AND.myid==0) write(io_lun,'(6x,i5,A)') myid, ' Zeroing matSFcoeff'
if(iprint_basis + min_layer>=5.AND.myid==0) write(io_lun,'(6x,i5,A)') myid, ' Zeroing matSFcoeff'
do spin_SF=1,nspin_SF
call matrix_scale(zero,matSFcoeff(spin_SF))
call matrix_scale(zero,matSFcoeff_tran(spin_SF))
enddo
if(iprint_basis>=5.AND.myid==0) write(io_lun,'(6x,A)') ' Done Zeroing'
if(iprint_basis + min_layer>=5.AND.myid==0) write(io_lun,'(6x,A)') ' Done Zeroing'
iprim = 0
call start_timer(tmr_std_matrices)
do part = 1,bundle%groups_on_node ! Loop over primary set partitions
if(iprint_basis>=6.AND.myid==0) write(io_lun,fmt='(6x,"Processor, partition: ",2i7)') myid,part
if(iprint_basis + min_layer>=6.AND.myid==0) write(io_lun,fmt='(6x,"Processor, partition: ",2i7)') myid,part
if(bundle%nm_nodgroup(part)>0) then ! If there are atoms in partition
do memb = 1,bundle%nm_nodgroup(part) ! Loop over atoms
atom_num = bundle%nm_nodbeg(part)+memb-1
iprim=iprim+1
! Atomic species
atom_spec = bundle%species(atom_num)
if(iprint_basis>=6.AND.myid==0) write(io_lun,'(6x,"Processor, atom, spec: ",3i7)') myid,memb,atom_spec
if(iprint_basis + min_layer>=6.AND.myid==0) write(io_lun,'(6x,"Processor, atom, spec: ",3i7)') myid,memb,atom_spec
do neigh = 1, mat(part,SFcoeff_range)%n_nab(memb) ! Loop over neighbours of atom
ist = mat(part,SFcoeff_range)%i_acc(memb)+neigh-1
! Build the distances between atoms - needed for phases
@ -427,7 +430,7 @@ contains
dz = BCS_parts%zcover(gcspart)-bundle%zprim(atom_num)
r2 = dx*dx + dy*dy + dz*dz
r2 = sqrt(r2)
if(iprint_basis>=6.AND.myid==0) write(80+myid,*) 'dx,y,z,r: ',gcspart,dx,dy,dz,r2
if(iprint_basis + min_layer>=6.AND.myid==0) write(80+myid,*) 'dx,y,z,r: ',gcspart,dx,dy,dz,r2
if (r2.le.RD_ERR) then ! only on-site
! We need to know the species of neighbour
neigh_global_part = BCS_parts%lab_cell(mat(part,SFcoeff_range)%i_part(ist))
@ -451,7 +454,7 @@ contains
enddo
! Write out current SF coefficients with some iprint (in future)
! if (iprint_basis>=3) call dump_pos_and_matrices
! if (iprint_basis + min_layer>=3) call dump_pos_and_matrices
! if (nspin_SF == 1) then
! call dump_matrix("SFcoeff", matSFcoeff(1), inode)
! else
@ -513,7 +516,7 @@ contains
! Loop over SFs on i
if (nsf_species(atom_spec).eq.npao_species(atom_spec)) then
! not-contracted SFs
if (iprint_basis>=6) write(io_lun,*) 'SFs of species', atom_spec, ' are not contracted'
if (iprint_basis + min_layer>=6) write(io_lun,*) 'SFs of species', atom_spec, ' are not contracted'
do sf1 = 1, nsf_i
count_pao_j = 1
! Loop over PAOs on j = i
@ -534,7 +537,7 @@ contains
enddo ! sf1
else
! contracted SFs
if (iprint_basis>=6) write(io_lun,*) 'SFs of species', atom_spec, ' are contracted'
if (iprint_basis + min_layer>=6) write(io_lun,*) 'SFs of species', atom_spec, ' are contracted'
do sf1 = 1, nsf_i
count_pao_j = 1 ! which PAO
! Loop over PAOs on j = i
@ -635,7 +638,7 @@ contains
real(double) :: r_center, r_shift, WIDTH
if (iprint_basis>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:smear_MSSF'
if (iprint_basis + min_layer>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:smear_MSSF'
sf1 = mat_p(matA)%sf1_type
sf2 = mat_p(matA)%sf2_type
@ -705,7 +708,7 @@ contains
else if (itype.eq.2) then
FLTR = half * erfc_cq(r2 - (r_center + r_shift)) ! f(r) = 0.5 * erfc(r - (r_center+r_shift))
endif
if (iprint_basis>=6) write(io_lun,'(2(A,F10.5))') 'r=',r2,' FLTR =',FLTR
if (iprint_basis + min_layer>=6) write(io_lun,'(2(A,F10.5))') 'r=',r2,' FLTR =',FLTR
do k = 1, nsf1
do l = 1, nsf2
@ -749,44 +752,48 @@ contains
use matrix_data, only: mat, SFcoeff_range
use global_module, only: ni_in_cell
use GenComms, only: gsum
use io_module, only: return_prefix
implicit none
! Local
integer :: stat, iprim, atom_i, np, i
integer, allocatable, dimension(:) :: num_naba_MS ! storage for number of naba atoms
character(len=14) :: subname = "Init_SFcoeff: "
character(len=120) :: prefix
if (iprint_basis + min_layer>2) then
prefix = return_prefix(subname, min_layer)
allocate(num_naba_MS(ni_in_cell),STAT=stat)
if(stat/=0) call cq_abort('print_naba_in_MSSF: error allocating num_naba_MS')
allocate(num_naba_MS(ni_in_cell),STAT=stat)
if(stat/=0) call cq_abort('print_naba_in_MSSF: error allocating num_naba_MS')
num_naba_MS(:) = 0
num_naba_MS(:) = 0
iprim = 0; atom_i = 0
do np = 1,bundle%groups_on_node ! Loop over primary set partitions
if(bundle%nm_nodgroup(np)>0) then ! If there are atoms in partition
do i = 1,bundle%nm_nodgroup(np) ! Loop over atom_i
iprim = iprim + 1
atom_i = bundle%ig_prim(iprim) ! global number of i
num_naba_MS(atom_i) = mat(np,SFcoeff_range)%n_nab(i)
enddo ! i
endif ! endif of (bundle%nm_nodgroup(np)>0)
enddo ! np
iprim = 0; atom_i = 0
do np = 1,bundle%groups_on_node ! Loop over primary set partitions
if(bundle%nm_nodgroup(np)>0) then ! If there are atoms in partition
do i = 1,bundle%nm_nodgroup(np) ! Loop over atom_i
iprim = iprim + 1
atom_i = bundle%ig_prim(iprim) ! global number of i
num_naba_MS(atom_i) = mat(np,SFcoeff_range)%n_nab(i)
enddo ! i
endif ! endif of (bundle%nm_nodgroup(np)>0)
enddo ! np
! print out the number of neighbour atoms for each target atom
call gsum(num_naba_MS,ni_in_cell)
if (inode==ionode .and. iprint_basis + min_layer>2) then
write(io_lun,fmt='(4x,a)') trim(prefix)//' --- Number of neighbour atoms in the multi-site range ---'
write(io_lun,fmt='(4x,a)') trim(prefix)//' atom ID number of neighbour atoms'
do i = 1, ni_in_cell
write(io_lun,'(4x,a,5X,I8,5X,I8)') trim(prefix),i, num_naba_MS(i)
enddo
endif
! print out the number of neighbour atoms for each target atom
call my_barrier()
call gsum(num_naba_MS,ni_in_cell)
if (inode==ionode .and. iprint_basis>1) then
write(io_lun,*) ' --- Number of neighbour atoms in the multi-site range ---'
write(io_lun,*) ' atom ID number of neighbour atoms'
do i = 1, ni_in_cell
write(io_lun,'(5X,I8,5X,I8)') i, num_naba_MS(i)
enddo
endif
! deallocate subspace matrices their labels
deallocate(num_naba_MS,STAT=stat)
if(stat/=0) call cq_abort('print_naba_in_MSSF: error deallocating num_naba_MS')
! deallocate subspace matrices their labels
deallocate(num_naba_MS,STAT=stat)
if(stat/=0) call cq_abort('print_naba_in_MSSF: error deallocating num_naba_MS')
end if
!
return
end subroutine print_naba_in_MSSF
@ -840,7 +847,7 @@ contains
matSatomf, matHatomf, matSFcoeff, aLa_aHa_aLHa, aLa_aSa_aLSa
use global_module, only: ni_in_cell, numprocs, nspin, nspin_SF, sf, atomf, flag_SpinDependentSF
use species_module, only: npao_species
use io_module, only: get_file_name
use io_module, only: get_file_name, return_prefix
use input_module, only: io_assign, io_close
use GenComms, only: mtime, gsum
use timer_stdclocks_module, only: start_timer,stop_timer,tmr_std_allocation
@ -873,10 +880,13 @@ contains
real(double), allocatable, dimension(:) :: l_k_r2
integer, allocatable, dimension(:) :: IWORK, IFAIL ! scratch space
real(double), allocatable, dimension(:) :: WORK, WORK2
character(len=7) :: subname = "DoLFD: "
character(len=120) :: prefix
prefix = return_prefix(subname, min_layer)
call my_barrier
if (inode==ionode .and. iprint_basis>2) write(io_lun,*) 'Do Localized Filter Diagonalization'
if (inode==ionode .and. iprint_basis + min_layer>2) write(io_lun,fmt='(4x,a)') &
trim(prefix)//' Entering'
do spin_SF = 1, nspin_SF
call matrix_scale(zero,matSFcoeff(spin_SF))
@ -889,7 +899,7 @@ contains
abstol = 1.0e-300_double
! open debug file for TVEC and subspace MOs
if (iprint_basis>=6) then
if (iprint_basis + min_layer>=6) then
call get_file_name('TVECr',numprocs,inode,filename11) ! Build a filename based on node number for TVEC
call io_assign(lun11) ! Open file
open(unit=lun11,file=filename11)
@ -905,7 +915,7 @@ contains
len_kj_sub = (nhalo_LFD*nhalo_LFD+nhalo_LFD)/2 ! max. number of lower-triangle of halo-pair (k,j)
len_Sub = (max_npao_LFD*max_npao_LFD+max_npao_LFD)/2 ! max. number of lower-triangle of subspace matrix elements
if (iprint_basis>=5 .and. inode==ionode) then
if (iprint_basis + min_layer>=5 .and. inode==ionode) then
if (flag_LFD_useChemPsub) then
write(io_lun,*) ' LFD: Chemical potential = ', 'will be determined later'
else
@ -946,7 +956,7 @@ contains
! --- Start Localized Filter Diagonalisation Method for each primary ATOM ---
!
if(myid==0) t0 = mtime()
if (iprint_basis>=5.and.inode==ionode) &
if (iprint_basis + min_layer>=5.and.inode==ionode) &
write(io_lun,'(/A,I2)') 'Start atom loop in the LFD method for spin',spin
iprim = 0; atom_i = 0
@ -970,7 +980,7 @@ contains
nd3 = LFDhalo%ndimj(k_in_halo) ! number of PAOs on atom_k
len_Sub_i = len_Sub_i + nd3 ! dimension of subspace for atom_i
enddo ! k
if (iprint_basis>=6) write(io_lun,'(A,I8,A,I2,A,I3,A,I4)') &
if (iprint_basis + min_layer>=6) write(io_lun,'(A,I8,A,I2,A,I3,A,I4)') &
' atom_i=',atom_i, ' : NTVEC=',NTVEC, &
' nab=',mat(np,LD_range)%n_nab(i), ' len_Sub_i = ',len_Sub_i
!
@ -1050,7 +1060,7 @@ contains
Ssub_i(:,:) = zero
call DCOPY(len_Sub_i*len_Sub_i,WORK2,1,Ssub_i,1) ! Ssub_i was overwritten, so reset to the original.
! debug: print out local MOs
if (iprint_basis>=6) call LFD_debug_matrix(lun12,1,np,atom_i,EVAL,EVEC,len_Sub_i,NUMEIG, &
if (iprint_basis + min_layer>=6) call LFD_debug_matrix(lun12,1,np,atom_i,EVAL,EVEC,len_Sub_i,NUMEIG, &
n_naba_i_d,l_k_g,l_k_r2,l_kpao)
!
! --- (4) filteration: k = Csub_i * f(e) * Csub_i**T * Ssub_i * TVEC
@ -1073,7 +1083,7 @@ contains
if (NEsub0.gt.NEsub) NEsub = NEsub + one
INEsub = idint(NEsub) ! number of occupied orbitals
ChemP = (EVAL(INEsub)+EVAL(INEsub+1)) / two ! ChemP is set to the average of subspace HOMO and LUMO
if (iprint_basis>=6) write(io_lun,'(A,i8,A,F10.5)') ' Atom',atom_i,': ChemP of this subspace =',ChemP
if (iprint_basis + min_layer>=6) write(io_lun,'(A,i8,A,F10.5)') ' Atom',atom_i,': ChemP of this subspace =',ChemP
endif
call LFD_filter(TVEC,EVAL,ChemP,kT1,NUMEIG,NTVEC)
!
@ -1088,8 +1098,9 @@ contains
WORK(:) = zero
call transpose_2Dmat(TVEC,WORK,len_Sub_i,NTVEC)
if (.not.flag_SpinDependentSF .and. spin.eq.2) then
if(inode==ionode .and. iprint_basis>1) &
write(io_lun,*) 'Take average of matSFcoeff(1) and matSFcoeff(2) into matSFcoeff(1)'
if(inode==ionode .and. iprint_basis + min_layer>2) &
write(io_lun,fmt='(4x,a)') &
trim(prefix)//'Take average of matSFcoeff(1) and matSFcoeff(2) into matSFcoeff(1)'
matSFcoeff_2 = allocate_temp_matrix(SFcoeff_range,0,sf,atomf)
call matrix_scale(zero,matSFcoeff_2)
call LFD_put_TVEC_to_SFcoeff(np,i,mat(:,SFcoeff_range),mat_p(matSFcoeff_2)%matrix, &
@ -1146,16 +1157,16 @@ contains
call stop_timer(tmr_std_allocation)
! Close debug file
if (iprint_basis>=6) then
if (iprint_basis + min_layer>=6) then
call io_close(lun11)
call io_close(lun12)
endif
if(myid==0 .and. iprint_basis>3) then
if(myid==0 .and. inode==ionode .and. iprint_basis + min_layer>2) then
t1 = mtime()
write(io_lun,'(A,f20.8,A)') 'Time for Localised Filter Diagonalisation: ',t1-t0,' ms'
write(io_lun,'(4x,A,f20.8,A)') trim(prefix)//' time for Localised Filter Diagonalisation: ',t1-t0,' ms'
end if
if (inode==ionode .and. iprint_basis>3) write(io_lun,*) 'Done Localized Filter Diagonalization'
if (inode==ionode .and. iprint_basis + min_layer>2) write(io_lun,fmt='(4x,a)') trim(prefix)//'Exiting'
!
return
end subroutine LocFilterDiag
@ -1253,7 +1264,7 @@ contains
integer :: count_kj
if (iprint_basis>=5 .and. inode==ionode) write(io_lun,*) 'We are in sub:LFD_make_Subspace_halo'
if (iprint_basis + min_layer>=5 .and. inode==ionode) write(io_lun,*) 'We are in sub:LFD_make_Subspace_halo'
call start_timer(tmr_std_allocation)
if(iprint_mat>3.AND.myid==0) t0 = mtime()
@ -1639,7 +1650,7 @@ contains
ist, gcspart, k_in_halo, j, ist2, gcspart2, j_in_halo, nd2, nd3, n2, n3
if (inode==ionode .and. iprint_basis>=5) write(io_lun,*) 'We are in sub:LFD_make_Subspace_i'
if (inode==ionode .and. iprint_basis + min_layer>4) write(io_lun,*) 'We are in sub:LFD_make_Subspace_i'
kpao0 = 0
! Loop over atom_k, neighbour of atom_i
@ -1842,7 +1853,7 @@ contains
! real(double) :: MSSF_nonminimal_offset
if (iprint_basis>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:LFD_make_TVEC'
if (iprint_basis + min_layer>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:LFD_make_TVEC'
npao_i = npao_species(atom_spec)
@ -2006,7 +2017,7 @@ contains
integer :: IVEC,JVEC
if (iprint_basis>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:LFD_filter'
if (iprint_basis + min_layer>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:LFD_filter'
do IVEC = 1,NUMEIG
EXFRM = ( EIG(IVEC) - CHEMP ) * kT1 ! (e_i - mu) / kT
@ -2050,7 +2061,7 @@ contains
integer :: i,j
if (iprint_basis>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:transpose_2Dmat'
if (iprint_basis + min_layer>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:transpose_2Dmat'
do i = 1,LDA
do j = 1,LDB
@ -2110,7 +2121,7 @@ contains
integer :: nb,ist_k,ist_l,nd1,nd2,bsize,SFCOEFFposn,LFDposn
if (iprint_basis>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:LFD_put_TVEC_to_SFcoeff'
if (iprint_basis + min_layer>=5.and.inode==ionode) write(io_lun,*) 'We are in sub:LFD_put_TVEC_to_SFcoeff'
ist_l = LFDmat(np)%i_acc(i) ! index of atom_l in neighbor labelling for LD_range
SFCOEFFposn = MSmat(np)%nd_offset+MSmat(np)%i_nd_acc(i) ! starting position of matrix element for (atom_i,atom_k) in SFCOEFF matrix

153
src/pao_minimisation.module.f90
View File

@ -129,7 +129,7 @@ contains
subroutine vary_pao(n_support_iterations, fixed_potential, vary_mu, &
n_cg_L_iterations, L_tolerance, sc_tolerance, &
energy_tolerance, total_energy_last, &
expected_reduction, level)
expected_reduction)
use datatypes
use logicals
@ -146,12 +146,13 @@ contains
ni_in_cell, &
flag_self_consistent, &
id_glob, numprocs, area_minE, &
nspin, spin_factor, nspin_SF, flag_diagonalisation
nspin, spin_factor, nspin_SF, &
flag_diagonalisation, min_layer
use group_module, only: parts
use H_matrix_module, only: get_H_matrix
use S_matrix_module, only: get_S_matrix
use store_matrix, only: dump_pos_and_matrices, unit_MSSF_save
! use io_module, only: dump_matrix
use io_module, only: return_prefix
use support_spec_format, only: TestBasisGrads, TestTot, &
TestBoth, TestS, TestH
use DMMin, only: FindMinDM
@ -166,14 +167,15 @@ contains
matdSFcoeff, matdSFcoeff_e, &
matrix_scale, matrix_transpose
use multisiteSF_module, only: normalise_SFcoeff, n_dumpSFcoeff
use units, only: en_conv, en_units, energy_units
use SelfCon, only: new_SC_potl
implicit none
! Shared variables
logical :: vary_mu, fixed_potential, convergence_flag
integer :: n_cg_L_iterations
integer :: n_support_iterations
integer :: level
real(double) :: expected_reduction
real(double) :: total_energy_last, energy_tolerance, L_tolerance, &
sc_tolerance
@ -197,8 +199,8 @@ contains
character(len=10) :: subname = "vary_pao: "
character(len=120) :: prefix
prefix = return_prefix(subname, level)
if(inode==ionode .and. iprint_minE + level >= 0) &
prefix = return_prefix(subname, min_layer)
if(inode==ionode .and. iprint_minE + min_layer >= 0) &
write(io_lun,fmt='(/4x,a)') trim(prefix)//" Starting PAO optimisation"
reset_L = .true.
@ -226,7 +228,7 @@ contains
con_tolerance = sc_tolerance
tolerance = L_tolerance
if (inode == ionode .and. iprint_minE + level >= 1) &
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, fmt='(4x,a,2e12.5)') trim(prefix)//' Tolerances: ', &
con_tolerance, tolerance
@ -460,7 +462,8 @@ contains
last_step = 1.0D10
! now loop over search directions
do n_iterations = 1, n_support_iterations
if (inode == ionode) write (io_lun, 7) n_iterations
if (inode == ionode .and. iprint_minE + min_layer >= 1) &
write (io_lun, fmt='(/4x,a,i3)') trim(prefix)//" PAO iteration ", n_iterations
do spin_SF = 1, nspin_SF
! We need the last search direction for CG manipulations
call copy(length, search_direction(:,spin_SF), 1, last_sd(:,spin_SF), 1)
@ -474,8 +477,8 @@ contains
mat_p(matdSFcoeff_e(spin_SF))%matrix, 1)
call gsum(dN_dot_de)
call gsum(dN_dot_dN)
if (inode == ionode) &
write (io_lun, *) 'dN.de, dN.dN ', &
if (inode == ionode .and. iprint_minE + min_layer >= 3) &
write (io_lun, fmt='(4x,a,2e12.5)') trim(prefix)//' dN.de, dN.dN ', &
dN_dot_de, dN_dot_dN
call axpy(length, - (dN_dot_de / dN_dot_dN), &
mat_p(matdSFcoeff_e(spin_SF))%matrix, 1, search_direction(:,spin_SF), 1)
@ -486,7 +489,6 @@ contains
gg(spin_SF) = dgg(spin_SF)
dgg(spin_SF) = dot(length, search_direction(:,spin_SF), 1, Psd(:,spin_SF), 1)
call gsum(dgg(spin_SF))
if (inode == ionode) write (io_lun, '(A,I2,A,F15.7)') 'dgg(',spin_SF,') is ', dgg(spin_SF)
if (gg(spin_SF) /= zero) then
gamma(spin_SF) = dgg(spin_SF) / gg(spin_SF)
else
@ -494,7 +496,9 @@ contains
end if
!gamma = zero
! if (mod(n_iterations, 5) == 0) gamma = zero
if (inode == ionode) write (io_lun, '(A,I2,A,F15.7)') 'Gamma(',spin_SF,') is ', gamma(spin_SF)
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, '(4x,A,I2,A,2F10.6)') trim(prefix)//' For spin(',spin_SF, &
') dgg and gamma are ', dgg(spin_SF), gamma(spin_SF)
! Construct the actual search direction
call copy(length, Psd(:,spin_SF), 1, search_direction(:,spin_SF), 1)
@ -507,8 +511,8 @@ contains
mat_p(matdSFcoeff_e(spin_SF))%matrix, 1)
call gsum(dN_dot_de)
call gsum(dN_dot_dN)
if (inode == ionode) &
write (io_lun, *) 'dN.de, dN.dN ', &
if (inode == ionode .and. iprint_minE + min_layer >=3) &
write (io_lun, fmt='(4x,a,2e12.5)') trim(prefix)//' dN.de, dN.dN ', &
dN_dot_de, dN_dot_dN
call axpy(length, - (dN_dot_de / dN_dot_dN), &
mat_p(matdSFcoeff_e(spin_SF))%matrix, 1, search_direction(:,spin_SF), 1)
@ -516,11 +520,11 @@ contains
! Check this !
sum_0(spin_SF) = dot(length, mat_p(matdSFcoeff(spin_SF))%matrix, 1, search_direction(:,spin_SF), 1)
call gsum(sum_0(spin_SF))
if (inode == ionode) write (io_lun, '(A,I2,A,F15.7)') 'sum_0(',spin_SF,') is ', sum_0(spin_SF)
if (inode == ionode .and. iprint_minE + min_layer >=3) &
write (io_lun, '(4x,A,I2,A,F15.7)') trim(prefix)//' sum_0(',spin_SF,') is ', &
sum_0(spin_SF)
call my_barrier()
! minimise the energy (approximately) in this direction.
if (inode == ionode) write (io_lun, *) 'Minimise'
call my_barrier()
tmp(spin_SF) = dot(length, search_direction(:,spin_SF), 1, mat_p(matdSFcoeff(spin_SF))%matrix, 1)
call gsum(tmp(spin_SF))
@ -535,25 +539,44 @@ contains
vary_mu, n_cg_L_iterations, tolerance, &
con_tolerance, total_energy_0, &
expected_reduction, last_step, tmp)
flag_vary_basis = .true.
if (inode == ionode) write (io_lun, *) 'Returned !'
! Normalise and writeout
!call normalise_SFcoeff
!do spin_SF = 1,nspin_SF
! call matrix_scale(zero,matSFcoeff_tran(spin_SF))
! call matrix_transpose(matSFcoeff(spin_SF), matSFcoeff_tran(spin_SF))
!enddo
! Write out current SF coefficients every n_dumpSFcoeff, if n_dumpSFcoeff > 0)
if (n_dumpSFcoeff > 0 .and. mod(n_iterations,n_dumpSFcoeff) == 1) then
call dump_pos_and_matrices(index = unit_MSSF_save)
endif
flag_vary_basis = .true.
! Find change in energy for convergence
diff = total_energy_last - total_energy_0
dE_PAO = diff
if (abs(diff / total_energy_0) <= energy_tolerance) then
if (inode == ionode) write (io_lun, 18) total_energy_0
if(inode==ionode .and. iprint_minE + min_layer>=0) &
write(io_lun, fmt='(4x,a,i3,2(a,f15.7,a2))') trim(prefix)//" Iter: ", &
n_iterations," E: ",total_energy_0*en_conv,en_units(energy_units),&
" dE: ",diff*en_conv,en_units(energy_units)
if (inode == ionode .and. iprint_minE + min_layer >=-1) &
write (io_lun, fmt='(/4x,a,f15.7,a2,a,i3,a)') &
trim(prefix)//" Minimisation converged to ",total_energy_0*en_conv, &
en_units(energy_units), " after ",n_iterations," iterations"
convergence_flag = .true.
total_energy_last = total_energy_0
deallocate(search_direction, last_sd, Psd, STAT=stat)
if (stat /= 0) call cq_abort("vary_pao: Error dealloc mem")
call reg_dealloc_mem(area_minE, 3*length*nspin_SF, type_dbl)
return
else
if(inode==ionode .and. iprint_minE + min_layer>=0) &
write(io_lun, fmt='(4x,a,i3,2(a,f15.7,a2))') trim(prefix)//" Iter: ", &
n_iterations," E: ",total_energy_0*en_conv,en_units(energy_units),&
" dE: ",diff*en_conv,en_units(energy_units)
end if
! We need to assemble the gradient
@ -566,8 +589,9 @@ contains
do spin_SF = 1, nspin_SF
summ = dot(length, mat_p(matdSFcoeff(spin_SF))%matrix, 1, mat_p(matdSFcoeff(spin_SF))%matrix, 1)
call gsum(summ)
if (inode == ionode) &
write (io_lun, *) 'Dot prod of gradient (',spin_SF,'): ', summ
if (inode == ionode .and. iprint_minE + min_layer >= 3) &
write (io_lun, fmt='(4x,a,i1,a,e12.5)') trim(prefix)//' Dot prod of gradient (',&
spin_SF,'): ', summ
enddo
total_energy_last = total_energy_0
end do ! n_iterations
@ -578,10 +602,6 @@ contains
return
7 format(/20x,'------------ PAO Variation #: ',i5,' ------------',/)
18 format(///20x,'The minimisation has converged to a total energy:', &
//20x,' Total energy = ',f15.7)
end subroutine vary_pao
!!***
@ -1225,13 +1245,15 @@ contains
use S_matrix_module, only: get_S_matrix
use GenComms, only: gsum, my_barrier, cq_abort, inode, &
ionode
! Check on whether we need K found from L or whether we have it exactly
use global_module, only: ni_in_cell, &
area_minE, nspin, nspin_SF, flag_diagonalisation
area_minE, nspin, nspin_SF, flag_diagonalisation, &
iprint_minE, min_layer
use primary_module, only: bundle
use memory_module, only: reg_alloc_mem, reg_dealloc_mem, &
type_dbl
use multisiteSF_module, only: normalise_SFcoeff
use units, only: en_conv, en_units, energy_units
use io_module, only: return_prefix
implicit none
@ -1257,10 +1279,14 @@ contains
logical :: done = .false. ! flag of line minimisation
real(double), save :: kmin_last = zero
real(double), save :: dE = zero ! Use this to guess initial step ?
if (inode == ionode) &
write (io_lun, *) 'On entry to pao line_min, dE is ', &
dE, total_energy_0
character(len=14) :: subname = "line_min_pao: "
character(len=120) :: prefix
min_layer = min_layer - 1
prefix = return_prefix(subname, min_layer)
if (inode == ionode .and. iprint_minE + min_layer >= 1) &
write (io_lun, fmt='(4x,a,f15.7,a2)') trim(prefix)//' On entry to pao line_min, dE is ', &
dE*en_conv, en_units(energy_units)
n_atoms = bundle%n_prim
length = mat_p(matSFcoeff(1))%length
@ -1268,7 +1294,8 @@ contains
tmp(spin_SF) = dot(length, search_direction(:,spin_SF), 1, search_direction(:,spin_SF), 1)
call gsum(tmp(spin_SF))
! search_direction(:,spin_SF) = search_direction(:,spin_SF) / sqrt(tmp(spin_SF))
if (inode == ionode) write (io_lun, *) 'Searchdir: ', tmp(spin_SF)
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, fmt='(4x,a,f15.7)') trim(prefix)//' Searchdir: ', tmp(spin_SF)
enddo
! if (nspin_SF == 1) then
@ -1330,19 +1357,21 @@ contains
! 3. Get a new self-consistent potential and Hamiltonian
! I've not put a call to get_H_matrix here because it's
! currently in new_SC_potl
min_layer = min_layer - 1
call new_SC_potl(.false., con_tolerance, reset_L, &
fixed_potential, vary_mu, n_cg_L_iterations, &
tolerance, e3)
if (inode == ionode) &
min_layer = min_layer + 1
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, &
fmt='(2x,"In pao_min, iter ",i3," &
fmt='(4x,a,i3," &
&step and energy are ",f15.10,f25.10)') &
iter, k3, e3
if (inode == ionode) &
trim(prefix)//" In pao_min, iter ",iter, k3, e3
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, &
fmt='(" iter=", i3," k1, k2, k3, &
&= ", 5f15.8)') &
iter, k1, k2, k3
fmt='(4x,a, i3," k1, k2, k3, &
&= ", 3f15.8)') &
trim(prefix)//" iter=",iter, k1, k2, k3
if (e3 < e2) then ! We're still going down hill
k1 = k2
e1 = e2
@ -1365,11 +1394,12 @@ contains
(k1 * k1 - k2 * k2) * (e1 - e3)) / &
((k1 - k3) * (e1 - e2) - (k1 - k2) * (e1 - e3)))
kmin_last = kmin
if (inode == ionode) &
write (io_lun, *) 'In pao_min, bracketed - min from extrap: ', &
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, fmt='(4x,a,4f15.7)') &
trim(prefix)//' In pao_min, bracketed - min from extrap: ', &
k1, k2, k3, kmin
if (inode == ionode) &
write (io_lun, *) 'In pao_min, bracketed - energies: ', &
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, fmt='(4x,a,3f15.7)') trim(prefix)//' In pao_min, bracketed - energies: ', &
e1, e2, e3
! Change blips: start from blip0
@ -1382,6 +1412,7 @@ contains
call matrix_scale(zero,matSFcoeff_tran(spin_SF))
call matrix_transpose(matSFcoeff(spin_SF), matSFcoeff_tran(spin_SF))
enddo
iter = iter + 1
! Find new self-consistent energy
! 1. Get new S matrix
@ -1395,8 +1426,10 @@ contains
! 3. Get a new self-consistent potential and Hamiltonian
! I've not put a call to get_H_matrix here because it's currently
! in new_SC_potl
min_layer = min_layer - 1
call new_SC_potl(.false., con_tolerance, reset_L, fixed_potential,&
vary_mu, n_cg_L_iterations, tolerance, energy_out)
min_layer = min_layer + 1
! If the interpolation failed, go back to the previous "minimum"
if (energy_out > e2) then
@ -1411,6 +1444,7 @@ contains
call matrix_scale(zero,matSFcoeff_tran(spin_SF))
call matrix_transpose(matSFcoeff(spin_SF), matSFcoeff_tran(spin_SF))
enddo
iter = iter + 1
! Find new self-consistent energy
! 1. Get new S matrix
@ -1422,23 +1456,28 @@ contains
end if
reset_L = .false. !true.
! 3. Get a new self-consistent potential and Hamiltonian
min_layer = min_layer - 1
call new_SC_potl(.false., con_tolerance, reset_L, &
fixed_potential, vary_mu, n_cg_L_iterations, &
tolerance, energy_out)
min_layer = min_layer + 1
end if
if (inode == ionode) &
write (io_lun, fmt='(2x,"In pao_min, at exit energy is ",f25.10)') &
energy_out
if (inode == ionode .and. iprint_minE + min_layer >= 0) &
write(io_lun, fmt='(4x,a,i3,a)') trim(prefix)//" Found minimum after ",iter," iterations"
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, fmt='(4x,a,f25.10)') &
trim(prefix)//" At exit energy is ",energy_out*en_conv,en_units(energy_units)
dE = total_energy_0 - energy_out
if (inode == ionode) &
write (io_lun, fmt='(2x,"On exit from pao line_min, dE is ",f25.10,2f25.10)') &
dE, total_energy_0, energy_out
if (inode == ionode .and. iprint_minE + min_layer >= 2) &
write (io_lun, fmt='(4x,a,f25.10,2f25.10,a2)') &
trim(prefix)//" On exit, dE is ",dE*en_conv, &
total_energy_0*en_conv, energy_out*en_conv,en_units(energy_units)
total_energy_0 = energy_out
deallocate(data_PAO0, STAT=stat)
if (stat /= 0) call cq_abort("line_minimise_pao: Error dealloc mem")
call reg_dealloc_mem(area_minE, length*nspin_SF, type_dbl)
min_layer = min_layer + 1
return
end subroutine line_minimise_pao
@ -1674,14 +1713,4 @@ contains
!sum = dot(nsf*mat(part,Srange)%n_nab(memb),data_K(nsf1,1:,point:),1,data_dHloc(npao1,1,point),1)
!gradient(npao1,nsf1,iprim) = gradient(npao1,nsf1,iprim) + sum
function return_prefix(name, level)
character(len=120) :: return_prefix
character(len=*) :: name
character(len=10):: prefix = " "
integer :: level
return_prefix = prefix(1:level)//name
end function return_prefix
end module pao_minimisation

23
src/pseudo_tm_module.f90
View File

@ -117,6 +117,7 @@ contains
subroutine init_pseudo_tm(ecore, ncf_out)
use datatypes, only: double
use numbers, only: zero
use global_module, only: iprint_pseudo, flag_neutral_atom
use block_module, only: n_pts_in_block
use GenComms, only: cq_abort, myid, inode, ionode
@ -175,7 +176,11 @@ contains
call stop_timer(tmr_std_pseudopot) ! This is timed inside set_tm_pseudo
call set_tm_pseudo !contained in this module
call start_timer(tmr_std_pseudopot)
call get_energy_shift(ecore) !contained in this module
if(flag_neutral_atom) then
ecore = zero
else
call get_energy_shift(ecore) !contained in this module
end if
call stop_timer(tmr_std_pseudopot)
!****lat<$
@ -1912,17 +1917,17 @@ contains
call calc_energy_shift(pseudo(ispecies),nfac2, eshift2, vlocG0_2)
call calc_energy_shift(pseudo(ispecies),nfac3, eshift3, vlocG0_3)
!if(abs(eshift2-eshift3) > eps) write(io_lun,fmt='(10x,"WARNING!!! eshift")')
if(inode == ionode.AND.iprint_pseudo>0) then
if(inode == ionode.AND.iprint_pseudo>2) then
write (io_lun, &
'(10x,a11,i2,a12,3i3,a19,d15.7,/,63x,d15.7,/,63x,d15.7)') &
'(4x,a11,i2,a12,3i3,a19,d15.7,/,63x,d15.7,/,63x,d15.7)') &
'ispecies = ', ispecies, ' eshift for', &
nfac1, nfac2, nfac3, ' times fine mesh = ', &
eshift1, eshift2, eshift3
write (io_lun,'(55x,a8,d15.7)') ' diff = ', eshift2 - eshift3
end if
if(inode == ionode.AND.iprint_pseudo>0) then
if(inode == ionode.AND.iprint_pseudo>2) then
write (io_lun, &
'(10x,a11,i2,a12,3i3,a19,d15.7,/,63x,d15.7,/,63x,d15.7)') &
'(4x,a11,i2,a12,3i3,a19,d15.7,/,63x,d15.7,/,63x,d15.7)') &
'ispecies = ', ispecies, ' eshift for', &
nfac1, nfac2, nfac3, ' times fine mesh = ', &
vlocG0_1, vlocG0_2, vlocG0_3
@ -2098,8 +2103,8 @@ contains
intvloc = intvloc + four *pi * rr * wos(imesh) * &
(vlocal + z * erfarg )
!(vlocal + z * erf(sqrt(beta) *rr) )
if(imesh==nmesh_vloc.AND.inode==ionode.AND.iprint_pseudo>1) &
write(io_lun,fmt='(2x,"Sum of local potential and core charge: ",3f12.8)') &
if(imesh==nmesh_vloc.AND.inode==ionode.AND.iprint_pseudo>2) &
write(io_lun,fmt='(4x,"Sum of local potential and core charge: ",3f12.8)') &
vlocal + z, vlocal, z
enddo
call start_timer(tmr_std_allocation)
@ -2109,7 +2114,7 @@ contains
!write(io_lun,*) 'Int. vloc: ',intvloc
!intvloc = intvloc + pi*z/beta
if(PRESENT(vlocG0)) vlocG0 = intvloc
!write(io_lun,*) 'Int. vloc, eshift: ',intvloc,eshift, beta, pi*z/beta, z/beta
!write(io_lun,*) 'Int. vloc, eshift: ',intvloc,eshift!, beta, pi*z/beta, z/beta
else
nmesh = pseudo%chlocal%n
nmesh_vloc = (nmesh-1) * nfac + 1
@ -2170,7 +2175,7 @@ contains
func2%f(imesh) = rho* four * pi * rr
zcore=zcore + func1%f(imesh)*wos(imesh)
enddo
if(inode == ionode.AND.iprint_pseudo>1) write(io_lun,fmt='(2x,"Integral of interpolated chlocal ",f8.3)') zcore
if(inode == ionode.AND.iprint_pseudo>2) write(io_lun,fmt='(4x,"Integral of interpolated chlocal ",f8.3)') zcore
zcore = - zcore ! makes zcore positive
!Makes vlocal & integrates vlocal