TABLE OF CONTENTS
- ABINIT/aprxdr
- ABINIT/dotprodm_v
- ABINIT/dotprodm_vn
- ABINIT/findminscf
- ABINIT/m_ab7_mixing
- ABINIT/scfeig
- ABINIT/sqnormm_v
- m_ab7_mixing/ab7_mixing_copy_current_step
- m_ab7_mixing/ab7_mixing_deallocate
- m_ab7_mixing/ab7_mixing_eval
- m_ab7_mixing/ab7_mixing_eval_allocate
- m_ab7_mixing/ab7_mixing_eval_deallocate
- m_ab7_mixing/ab7_mixing_new
- m_ab7_mixing/ab7_mixing_use_disk_cache
- m_ab7_mixing/ab7_mixing_use_moving_atoms
- m_ab7_mixing/init_
- m_ab7_mixing/nullify
- m_ab7_mixing/scfcge
- m_ab7_mixing/scfopt
ABINIT/aprxdr [ Functions ]
NAME
aprxdr
FUNCTION
Compute the approximative derivatives of the energy at different points along the line search, thanks to a finite-difference formula. This formula is the projection along the line search of the Eq.(11) in PRB54, 4383 (1996) [[cite:Gonze1996]].
INPUTS
cplex: if 1, real space functions on FFT grid are REAL, if 2, COMPLEX choice= if==3, compute dedv_new, dedv_old, and dedv_mix, if/=3, compute only dedv_new and dedv_old. i_vresid and i_rhor, see the next lines. f_fftgr(nfft,nspden,n_fftgr)=different functions defined on the fft grid : The last residual potential is in f_fftgr(:,:,i_vresid(1)). The old residual potential is in f_fftgr(:,:,i_vresid(2)). The previous old residual potential is in f_fftgr(:,:,i_vresid(3)). (needed only when choice==3) The old density is in f_fftgr(:,:,i_rhor2). f_atm(3,natom,n_fftgr)=different functions defined for each atom : The last HF force is in f_atm(:,:,i_vresid(1)). The old HF force is in f_fftgr(:,:,i_vresid(2)). The previous old HF force is in f_fftgr(:,:,i_vresid(3)). (needed only when choice==3) The old atomic positions are in f_atm(:,:,i_rhor2) moved_atm_inside: if==1, the atoms are allowed to move. mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT natom=number of atoms in unit cell nfft=(effective) number of FFT grid points (for this processor) nfftot=total number of FFT grid points nspden=number of spin-density components rhor(nfft,nspden)=actual density xred(3,natom)=reduced atomic coordinates
OUTPUT
dedv_mix=approximate derivative from previous old residual dedv_new=approximate derivative from new residual dedv_old=approximate derivative from old residual (output only when choice==3)
NOTES
Should be OpenMP parallelized
SOURCE
2948 subroutine aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,& 2949 & f_atm,f_fftgr,i_rhor2,i_vresid,moved_atm_inside,& 2950 & mpicomm,mpi_summarize,natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred) 2951 2952 !Arguments ------------------------------------ 2953 !scalars 2954 integer,intent(in) :: choice,cplex,i_rhor2,moved_atm_inside,n_fftgr,natom,nfft 2955 integer,intent(in) :: mpicomm,nfftot,nspden 2956 logical, intent(in) :: mpi_summarize 2957 real(dp),intent(in) :: ucvol 2958 real(dp),intent(out) :: dedv_mix,dedv_new,dedv_old 2959 !arrays 2960 integer,intent(in) :: i_vresid(3) 2961 real(dp),intent(in) :: f_atm(3,natom,n_fftgr) 2962 real(dp),intent(in) :: f_fftgr(cplex*nfft,nspden,n_fftgr) 2963 real(dp),intent(in) :: rhor(cplex*nfft,nspden),xred(3,natom) 2964 2965 !Local variables------------------------------- 2966 !scalars 2967 integer :: iatom,idir 2968 !arrays 2969 real(dp) :: dedv_temp(1) 2970 real(dp),allocatable :: ddens(:,:,:) 2971 2972 ! ************************************************************************* 2973 2974 ABI_MALLOC(ddens,(cplex*nfft,nspden,1)) 2975 2976 !Compute approximative derivative of the energy 2977 !with respect to change of potential 2978 2979 ddens(:,:,1)=rhor(:,:)-f_fftgr(:,:,i_rhor2) 2980 2981 !call dotprod_vn(cplex,1,ddens,dedv_old,nfft,nfftot,nspden,1,vresid,ucvol) 2982 !Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(2)) 2983 call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(2),mpicomm,mpi_summarize,1,1,1,& 2984 & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol) 2985 dedv_old = dedv_temp(1) 2986 2987 !Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(1)) 2988 call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(1),mpicomm,mpi_summarize,1,1,1,& 2989 & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol) 2990 dedv_new= dedv_temp(1) 2991 2992 if(choice==3)then 2993 ! Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(3)) 2994 call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(3),mpicomm,mpi_summarize,1,1,1,& 2995 & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol) 2996 dedv_mix = dedv_temp(1) 2997 end if 2998 2999 ABI_FREE(ddens) 3000 3001 !------------------------------------------------------- 3002 3003 !Now, take care of eventual atomic displacements 3004 3005 if(moved_atm_inside==1)then 3006 do idir=1,3 3007 do iatom=1,natom 3008 dedv_new=dedv_new+& 3009 & f_atm(idir,iatom,i_vresid(1))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2)) 3010 dedv_old=dedv_old+& 3011 & f_atm(idir,iatom,i_vresid(2))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2)) 3012 if(choice==3) dedv_mix=dedv_mix+& 3013 & f_atm(idir,iatom,i_vresid(3))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2)) 3014 end do 3015 end do 3016 end if 3017 3018 end subroutine aprxdr
ABINIT/dotprodm_v [ Functions ]
NAME
dotprodm_v
FUNCTION
For two sets of potentials, compute dot product of each pair of two potentials (integral over FFT grid), to obtain a series of square residual-like quantity (so the sum of product of values is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume). Take into account the spin components of the potentials (nspden), and sum over them. Need the index of the first pair of potentials to be treated, in each array of potentials, and the number of potentials to be treated. Might be used to compute just one square of norm, in a big array, such as to avoid copying a potential from a big array to a temporary place.
INPUTS
cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX cpldot=if 1, the dot array is real, if 2, the dot array is complex index1=index of the first potential to be treated in the potarr1 array index2=index of the first potential to be treated in the potarr2 array mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT mult1=number of potentials to be treated in the first set mult2=number of potentials to be treated in the second set nfft= (effective) number of FFT grid points (for this processor) npot1= third dimension of the potarr1 array npot2= third dimension of the potarr2 array nspden=number of spin-density components opt_storage: 0, if potentials are stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn] 1, if potentials are stored as V, B_x, B_y, Bz (B=magn. field) potarr1(cplex*nfft,nspden,npot)=first array of real space potentials on FFT grid (if cplex=2 and cpldot=2, potarr1 is the array that will be complex conjugated) potarr2(cplex*nfft,nspden,npot)=second array of real space potentials on FFT grid
OUTPUT
dot(cpldot,mult1,mult2)= series of values of the dot product
SIDE EFFECTS
NOTES
Concerning storage when nspden=4:
cplex=1:
opt_storage=0: V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
opt_storage=1: V are stored as : V, B_x, B_y, B_z (real)
cplex=2:
opt_storage=0: V are stored as : V^11, V^22, V^12, i.V^21 (complex)
opt_storage=1: V are stored as : V, B_x, B_y, B_z (complex)
SOURCE
2382 subroutine dotprodm_v(cplex,cpldot,dot,index1,index2,mpicomm,mpi_summarize,& 2383 & mult1,mult2,nfft,npot1,npot2,nspden,opt_storage,potarr1,potarr2) 2384 2385 !Arguments ------------------------------------ 2386 !scalars 2387 integer,intent(in) :: cpldot,cplex,index1,index2,mult1,mult2,nfft,npot1,npot2 2388 integer,intent(in) :: nspden,opt_storage,mpicomm 2389 logical, intent(in) :: mpi_summarize 2390 !arrays 2391 real(dp),intent(in) :: potarr1(cplex*nfft,nspden,npot1) 2392 real(dp),intent(in) :: potarr2(cplex*nfft,nspden,npot2) 2393 real(dp),intent(out) :: dot(cpldot,mult1,mult2) 2394 2395 !Local variables------------------------------- 2396 !scalars 2397 integer :: i1,i2,ierr,ifft,ispden 2398 real(dp) :: ai,ar 2399 !arrays 2400 real(dp) :: tsec(2) 2401 2402 ! ************************************************************************* 2403 2404 !Real or complex inputs are coded 2405 DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex") 2406 2407 !Real or complex outputs are coded 2408 DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot") 2409 DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden") 2410 DBG_CHECK( npot1-index1-mult1 >= -1,"npot1-index1-mult1") 2411 DBG_CHECK( npot2-index2-mult2 >= -1,"npot2-index2-mult2") 2412 2413 if(cplex==1 .or. cpldot==1)then 2414 2415 do i1=1,mult1 2416 do i2=1,mult2 2417 ar=zero 2418 do ispden=1,min(nspden,2) 2419 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar) 2420 do ifft=1,cplex*nfft 2421 ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1) 2422 end do 2423 end do 2424 dot(1,i1,i2)=ar 2425 if (nspden==4) then 2426 ar=zero 2427 do ispden=3,4 2428 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar) 2429 do ifft=1,cplex*nfft 2430 ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1) 2431 end do 2432 end do 2433 if (opt_storage==0) then 2434 if (cplex==1) then 2435 dot(1,i1,i2)=dot(1,i1,i2)+two*ar 2436 else 2437 dot(1,i1,i2)=dot(1,i1,i2)+ar 2438 end if 2439 else 2440 dot(1,i1,i2)=half*(dot(1,i1,i2)+ar) 2441 end if 2442 end if 2443 end do 2444 end do 2445 2446 else ! if (cplex==2 .and. cpldot==2) 2447 2448 do i1=1,mult1 2449 do i2=1,mult2 2450 ar=zero ; ai=zero 2451 do ispden=1,min(nspden,2) 2452 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai) 2453 do ifft=1,nfft 2454 ar=ar + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) & 2455 & + potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1) 2456 ai=ai + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1) & 2457 & - potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) 2458 end do 2459 end do 2460 dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai 2461 if (nspden==4) then 2462 ar=zero 2463 do ispden=3,4 2464 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai) 2465 do ifft=1,nfft 2466 ar=ar + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) & 2467 & + potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1) 2468 ai=ai + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1) & 2469 & - potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) 2470 end do 2471 end do 2472 if (opt_storage==0) then 2473 dot(1,i1,i2)=dot(1,i1,i2)+ar 2474 dot(2,i1,i2)=dot(2,i1,i2)+ai 2475 else 2476 dot(1,i1,i2)=half*(dot(1,i1,i2)+ar) 2477 dot(2,i1,i2)=half*(dot(2,i1,i2)+ai) 2478 end if 2479 end if 2480 end do 2481 end do 2482 end if 2483 2484 !XG030513 : MPIWF reduction (addition) on dot is needed here 2485 if (mpi_summarize) then 2486 call timab(48,1,tsec) 2487 call xmpi_sum(dot,mpicomm ,ierr) 2488 call timab(48,2,tsec) 2489 end if 2490 2491 if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero 2492 2493 end subroutine dotprodm_v
ABINIT/dotprodm_vn [ Functions ]
NAME
dotprodm_vn
FUNCTION
For a set of densities and a set of potentials, compute the dot product (integral over FFT grid) of each pair, to obtain a series of energy-like quantity (so the usual dotproduct is divided by the number of FFT points, and multiplied by the primitive cell volume). Take into account the spin components of the density and potentials (nspden), and sum correctly over them. Note that the storage of densities and potentials is different : for potential, one stores the matrix components, while for the density, one stores the trace, and then, either the spin-polarisation (if nspden=2), or the magnetization vector (if nspden=4). Need the index of the first density/potential pair to be treated, in each array, and the number of pairs to be treated. Might be used to compute just one dot product, in a big array, such as to avoid copying the density and potential from a big array to a temporary place.
INPUTS
cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX cpldot=if 1, the dot array is real, if 2, the dot array is complex (not coded yet for nspden=4) denarr(cplex*nfft,nspden,nden)=real space density on FFT grid id=index of the first density to be treated in the denarr array ip=index of the first potential to be treated in the potarr array mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT multd=number of densities to be treated multp=number of potentials to be treated nden=third dimension of the denarr array nfft= (effective) number of FFT grid points (for this processor) nfftot= total number of FFT grid points npot=third dimension of the potarr array nspden=number of spin-density components potarr(cplex*nfft,nspden,npot)=real space potential on FFT grid (will be complex conjugated if cplex=2 and cpldot=2) ucvol=unit cell volume (Bohr**3)
OUTPUT
dot(cpldot,multp,multd)= series of values of the dot product potential/density
SIDE EFFECTS
NOTES
Concerning storage when nspden=4:
cplex=1:
V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
N are stored as : n, m_x, m_y, m_z (real)
cplex=2:
V are stored as : V^11, V^22, V^12, i.V^21 (complex)
N are stored as : n, m_x, m_y, mZ (complex)
SOURCE
2551 subroutine dotprodm_vn(cplex,cpldot,denarr,dot,id,ip,mpicomm, mpi_summarize,multd,multp,& 2552 & nden,nfft,nfftot,npot,nspden,potarr,ucvol) 2553 2554 !Arguments ------------------------------------ 2555 !scalars 2556 integer,intent(in) :: cpldot,cplex,id,ip,multd,multp,nden,nfft,nfftot,npot 2557 integer,intent(in) :: nspden,mpicomm 2558 logical, intent(in) :: mpi_summarize 2559 real(dp),intent(in) :: ucvol 2560 !arrays 2561 real(dp),intent(in) :: denarr(cplex*nfft,nspden,nden) 2562 real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot) 2563 real(dp),intent(out) :: dot(cpldot,multp,multd) 2564 2565 !Local variables------------------------------- 2566 !scalars 2567 integer :: i1,i2,ierr,ir,jr 2568 real(dp) :: ai,ar,dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21 2569 real(dp) :: dre22,dre_dn,dre_up,factor,pim11,pim12,pim21,pim22,pim_dn,pim_up 2570 real(dp) :: pre11,pre12,pre21,pre22,pre_dn,pre_up 2571 !arrays 2572 real(dp) :: tsec(2) 2573 2574 ! ************************************************************************* 2575 2576 !Real or complex inputs are coded 2577 DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex") 2578 2579 !Real or complex outputs are coded 2580 DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot") 2581 DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden") 2582 2583 DBG_CHECK(id >= 1,'Wrong id') 2584 DBG_CHECK(ip >= 1,'Wrong id') 2585 2586 DBG_CHECK(multd >= 1,"wrong multd") 2587 DBG_CHECK(multp >= 1,"wrong multp") 2588 2589 DBG_CHECK(nden-id-multd >=-1,'nden-id-multd') 2590 DBG_CHECK(npot-ip-multp >=-1,'npot-ip-multp') 2591 2592 if(nspden==1)then 2593 2594 if(cpldot==1 .or. cplex==1 )then 2595 2596 do i2=1,multd 2597 do i1=1,multp 2598 ar=zero 2599 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar) 2600 do ir=1,cplex*nfft 2601 ar=ar + potarr(ir,1,ip+i1-1)*denarr(ir,1,id+i2-1) 2602 end do 2603 dot(1,i1,i2)=ar 2604 end do ! i1 2605 end do ! i2 2606 2607 else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot 2608 2609 do i2=1,multd 2610 do i1=1,multp 2611 ar=zero ; ai=zero 2612 !$OMP PARALLEL DO PRIVATE(ir,jr) SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai) 2613 do ir=1,nfft 2614 jr=2*ir 2615 ar=ar + potarr(jr-1,1,ip+i1-1)*denarr(jr-1,1,id+i2-1) & 2616 & + potarr(jr ,1,ip+i1-1)*denarr(jr ,1,id+i2-1) 2617 ai=ai + potarr(jr-1,1,ip+i1-1)*denarr(jr ,1,id+i2-1) & 2618 & - potarr(jr ,1,ip+i1-1)*denarr(jr-1,1,id+i2-1) 2619 end do 2620 dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai 2621 end do ! i1 2622 end do ! i2 2623 2624 end if 2625 2626 else if(nspden==2)then 2627 2628 if(cpldot==1 .or. cplex==1 )then 2629 2630 do i2=1,multd 2631 do i1=1,multp 2632 ar=zero 2633 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar) 2634 do ir=1,cplex*nfft 2635 ar=ar + potarr(ir,1,ip+i1-1)* denarr(ir,2,id+i2-1) & ! This is the spin up contribution 2636 & + potarr(ir,2,ip+i1-1)*(denarr(ir,1,id+i2-1)-denarr(ir,2,id+i2-1)) ! This is the spin down contribution 2637 end do 2638 dot(1,i1,i2)=ar 2639 end do ! i1 2640 end do ! i2 2641 2642 else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot 2643 2644 do i2=1,multd 2645 do i1=1,multp 2646 ar=zero ; ai=zero 2647 !$OMP PARALLEL DO PRIVATE(ir,jr,dre_up,dim_up,dre_dn,dim_dn,pre_up,pim_up,pre_dn,pim_dn) & 2648 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai) 2649 do ir=1,nfft 2650 jr=2*ir 2651 2652 dre_up=denarr(jr-1,2,id+i2-1) 2653 dim_up=denarr(jr ,2,id+i2-1) 2654 dre_dn=denarr(jr-1,1,id+i2-1)-dre_up 2655 dim_dn=denarr(jr ,1,id+i2-1)-dim_up 2656 2657 pre_up=potarr(jr-1,1,ip+i1-1) 2658 pim_up=potarr(jr ,1,ip+i1-1) 2659 pre_dn=potarr(jr-1,2,ip+i1-1) 2660 pim_dn=potarr(jr ,2,ip+i1-1) 2661 2662 ar=ar + pre_up * dre_up & 2663 & + pim_up * dim_up & 2664 & + pre_dn * dre_dn & 2665 & + pim_dn * dim_dn 2666 ai=ai + pre_up * dim_up & 2667 & - pim_up * dre_up & 2668 & + pre_dn * dim_dn & 2669 & - pim_dn * dre_dn 2670 2671 end do 2672 dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai 2673 end do ! i1 2674 end do ! i2 2675 2676 end if 2677 2678 else if(nspden==4)then 2679 ! \rho{\alpha,\beta} V^{\alpha,\beta} = 2680 ! rho*(V^{11}+V^{22})/2$ 2681 ! + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2 2682 if (cplex==1) then 2683 do i2=1,multd 2684 do i1=1,multp 2685 ar=zero 2686 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar) 2687 do ir=1,cplex*nfft 2688 ar=ar+(potarr(ir,1,ip+i1-1)+potarr(ir,2,ip+i1-1))*half*denarr(ir,1,id+i2-1)& ! This is the density contrib 2689 & + potarr(ir,3,ip+i1-1) *denarr(ir,2,id+i2-1)& ! This is the m_x contrib 2690 & - potarr(ir,4,ip+i1-1) *denarr(ir,3,id+i2-1)& ! This is the m_y contrib 2691 & +(potarr(ir,1,ip+i1-1)-potarr(ir,2,ip+i1-1))*half*denarr(ir,4,id+i2-1) ! This is the m_z contrib 2692 end do 2693 dot(1,i1,i2)=ar 2694 end do ! i1 2695 end do ! i2 2696 else ! cplex=2 2697 ! Note concerning storage when cplex=2: 2698 ! V are stored as : v^11, v^22, V^12, i.V^21 (each are complex) 2699 ! N are stored as : n, m_x, m_y, mZ (each are complex) 2700 if (cpldot==1) then 2701 do i2=1,multd 2702 do i1=1,multp 2703 ar=zero ; ai=zero 2704 !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre22,pim22,pre12,pim12) & 2705 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar) 2706 do ir=1,nfft 2707 jr=2*ir 2708 dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2)) 2709 dim11=half*(denarr(jr ,1,id+i2)+denarr(jr-1,4,id+i2)) 2710 dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2)) 2711 dim22=half*(denarr(jr ,1,id+i2)-denarr(jr-1,4,id+i2)) 2712 dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr ,3,id+i2)) 2713 dim12=half*(denarr(jr ,2,id+i2)-denarr(jr-1,3,id+i2)) 2714 dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr ,3,id+i2)) 2715 dim21=half*(denarr(jr ,2,id+i2)+denarr(jr-1,3,id+i2)) 2716 pre11= potarr(jr-1,1,ip+i1) 2717 pim11= potarr(jr ,1,ip+i1) 2718 pre22= potarr(jr-1,2,ip+i1) 2719 pim22= potarr(jr ,2,ip+i1) 2720 pre12= potarr(jr-1,3,ip+i1) 2721 pim12= potarr(jr ,3,ip+i1) 2722 pre21= potarr(jr ,4,ip+i1) 2723 pim21=-potarr(jr-1,4,ip+i1) 2724 ar=ar + pre11 * dre11 & 2725 & + pim11 * dim11 & 2726 & + pre22 * dre22 & 2727 & + pim22 * dim22 & 2728 & + pre12 * dre12 & 2729 & + pim12 * dim12 & 2730 & + pre21 * dre21 & 2731 & + pim21 * dim21 2732 end do 2733 dot(1,i1,i2)=ar 2734 end do ! i1 2735 end do ! i2 2736 else !cpldot=2 2737 do i2=1,multd 2738 do i1=1,multp 2739 ar=zero ; ai=zero 2740 !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre12,pim12,pre22,pim22) & 2741 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai) 2742 do ir=1,nfft 2743 jr=2*ir 2744 dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2)) 2745 dim11=half*(denarr(jr ,1,id+i2)+denarr(jr-1,4,id+i2)) 2746 dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2)) 2747 dim22=half*(denarr(jr ,1,id+i2)-denarr(jr-1,4,id+i2)) 2748 dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr ,3,id+i2)) 2749 dim12=half*(denarr(jr ,2,id+i2)-denarr(jr-1,3,id+i2)) 2750 dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr ,3,id+i2)) 2751 dim21=half*(denarr(jr ,2,id+i2)+denarr(jr-1,3,id+i2)) 2752 pre11= potarr(jr-1,1,ip+i1) 2753 pim11= potarr(jr ,1,ip+i1) 2754 pre22= potarr(jr-1,2,ip+i1) 2755 pim22= potarr(jr ,2,ip+i1) 2756 pre12= potarr(jr-1,3,ip+i1) 2757 pim12= potarr(jr ,3,ip+i1) 2758 pre21= potarr(jr ,4,ip+i1) 2759 pim21=-potarr(jr-1,4,ip+i1) 2760 ar=ar + pre11 * dre11 & 2761 & + pim11 * dim11 & 2762 & + pre22 * dre22 & 2763 & + pim22 * dim22 & 2764 & + pre12 * dre12 & 2765 & + pim12 * dim12 & 2766 & + pre21 * dre21 & 2767 & + pim21 * dim21 2768 ai=ai + pre11 * dim11 & 2769 & - pim11 * dre11 & 2770 & + pre22 * dim22 & 2771 & - pim22 * dre22 & 2772 & + pre12 * dim12 & 2773 & - pim12 * dre12 & 2774 & + pre21 * dim21 & 2775 & - pim21 * dre21 2776 end do 2777 dot(1,i1,i2)=ar 2778 dot(2,i1,i2)=ai 2779 end do ! i1 2780 end do ! i2 2781 end if ! cpldot 2782 end if ! cplex 2783 end if ! nspden 2784 2785 factor=ucvol/dble(nfftot) 2786 dot(:,:,:)=factor*dot(:,:,:) 2787 2788 !XG030513 : MPIWF reduction (addition) on dot is needed here 2789 if (mpi_summarize) then 2790 call timab(48,1,tsec) 2791 call xmpi_sum(dot,mpicomm ,ierr) 2792 call timab(48,2,tsec) 2793 end if 2794 2795 if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero 2796 2797 end subroutine dotprodm_vn
ABINIT/findminscf [ Functions ]
NAME
findminscf
FUNCTION
Compute the minimum of a function whose value and derivative are known at two points, using different algorithms. Also deduce different quantities at this predicted point, and at the two other points
INPUTS
choice=1,uses a linear interpolation of the derivatives
=2,uses a quadratic interpolation based on the
values of the function, and the second derivative at mid-point
etotal_1=first value of the function
etotal_2=second value of the function
dedv_1=first value of the derivative
dedv_2=second value of the derivative
lambda_1=first value of the argument
lambda_2=second value of the argument
OUTPUT
dedv_predict=predicted value of the derivative (usually zero,
except if choice=4, if it happens that a minimum cannot be located,
and a trial step is taken)
d2edv2_predict=predicted value of the second derivative (not if choice=4)
d2edv2_1=first value of the second derivative (not if choice=4)
d2edv2_2=second value of the second derivative (not if choice=4)
etotal_predict=predicted value of the function
lambda_predict=predicted value of the argument
status= 0 if everything went normally ;
1 if negative second derivative
2 if some other problem
SOURCE
2243 subroutine findminscf(choice,dedv_1,dedv_2,dedv_predict,& 2244 & d2edv2_1,d2edv2_2,d2edv2_predict,& 2245 & etotal_1,etotal_2,etotal_predict,& 2246 & lambda_1,lambda_2,lambda_predict,errid,errmess) 2247 2248 !Arguments ------------------------------------ 2249 !scalars 2250 integer,intent(in) :: choice 2251 integer,intent(out) :: errid 2252 character(len=500), intent(out) :: errmess 2253 real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2 2254 real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict 2255 real(dp),intent(out) :: etotal_predict,lambda_predict 2256 2257 !Local variables------------------------------- 2258 !scalars 2259 real(dp) :: cc,d2edv2_mid,d_lambda,dedv_2bis 2260 real(dp) :: dedv_mid2,etotal_2bis 2261 character(len=500) :: message 2262 2263 ! ************************************************************************* 2264 2265 !DEBUG 2266 !write(std_out,*)' findmin : enter' 2267 !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2 2268 !ENDDEBUG 2269 2270 errid = AB7_NO_ERROR 2271 d_lambda=lambda_1-lambda_2 2272 2273 if(choice==1) then 2274 2275 ! Use the derivative information to predict lambda 2276 d2edv2_mid=(dedv_1-dedv_2)/d_lambda 2277 lambda_predict=lambda_2-dedv_2/d2edv2_mid 2278 dedv_predict=dedv_2+(lambda_predict-lambda_2)*d2edv2_mid 2279 d2edv2_1=d2edv2_mid 2280 d2edv2_2=d2edv2_mid 2281 d2edv2_predict=d2edv2_mid 2282 ! also use the first energy to predict new energy 2283 etotal_predict=etotal_1+dedv_1*(lambda_predict-lambda_1)& 2284 & +0.5_dp*d2edv2_1*(lambda_predict-lambda_1)**2 2285 etotal_2bis=etotal_1+dedv_1*(lambda_2-lambda_1)& 2286 & +0.5_dp*d2edv2_1*(lambda_2-lambda_1)**2 2287 2288 if(d2edv2_mid<0.0_dp)then 2289 errid = AB7_ERROR_MIXING_INTERNAL 2290 write(errmess,'(a,es18.10,a)')'The second derivative is negative, equal to ',d2edv2_mid,'.' 2291 ABI_WARNING(errmess) 2292 end if 2293 2294 else if(choice==2) then 2295 2296 ! Use energies and first derivative information 2297 ! etotal = aa + bb * lambda + cc * lambda**2 2298 dedv_mid2=(etotal_1-etotal_2)/d_lambda 2299 cc=(dedv_1-dedv_mid2)/d_lambda 2300 lambda_predict=lambda_1-0.5_dp*dedv_1/cc 2301 d2edv2_1=2*cc 2302 d2edv2_2=d2edv2_1 2303 d2edv2_predict=d2edv2_1 2304 if(d2edv2_predict<0.0_dp)then 2305 errid = AB7_ERROR_MIXING_INTERNAL 2306 write(errmess, '(a,es18.10,a,a,a)' )& 2307 & 'The second derivative is negative, equal to',d2edv2_predict,'.',ch10,& 2308 & '=> Pivoting ' 2309 ABI_WARNING(errmess) 2310 if(etotal_2 < etotal_1)then 2311 lambda_predict=lambda_2-0.5_dp*(lambda_1-lambda_2) 2312 else 2313 lambda_predict=lambda_1-0.5_dp*(lambda_2-lambda_1) 2314 end if 2315 end if 2316 dedv_predict=dedv_1+(lambda_predict-lambda_1)*d2edv2_1 2317 dedv_2bis=dedv_1+(lambda_2-lambda_1)*d2edv2_1 2318 etotal_predict=etotal_1+dedv_1*(lambda_predict-lambda_1)& 2319 & +0.5_dp*d2edv2_1*(lambda_predict-lambda_1)**2 2320 2321 end if 2322 2323 write(message, '(a,es12.4,a,es18.10)' ) & 2324 & ' findmin : lambda_predict ',lambda_predict,' etotal_predict ',etotal_predict 2325 call wrtout(std_out,message,'COLL') 2326 2327 end subroutine findminscf
ABINIT/m_ab7_mixing [ Modules ]
NAME
m_ab7_mixing
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (XG, DC, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 22 module m_ab7_mixing 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_linalg_interfaces 28 use m_xmpi 29 30 use m_time, only : timab 31 use m_io_tools, only : open_file 32 33 implicit none 34 35 private
ABINIT/scfeig [ Functions ]
NAME
scfeig
FUNCTION
Compute the largest eigenvalue and eigenvector of the SCF cycle. A brute force algorithm is presently used.
INPUTS
istep= number of the step in the SCF cycle nfft=(effective) number of FFT grid points (for this processor) nspden=number of spin-density components
OUTPUT
(see side effects)
SIDE EFFECTS
vtrial0(nfft,nspden)= contains vtrial at istep == 1
vtrial(nfft,nspden)= at input, it is the trial potential that gave vresid .
at output, it is an updated trial potential
vrespc(nfft,nspden)=the input preconditioned residual potential
work(nfft,nspden,2)=work space
SOURCE
1683 subroutine scfeig(istep,nfft,nspden,vrespc,vtrial,vtrial0,work,errid,errmess) 1684 1685 !Arguments ------------------------------------ 1686 !scalars 1687 integer,intent(in) :: istep,nfft,nspden 1688 integer,intent(out) :: errid 1689 character(len = 500), intent(out) :: errmess 1690 !arrays 1691 real(dp),intent(inout) :: vtrial0(nfft,nspden),work(nfft,nspden,2) 1692 real(dp),intent(inout) :: vrespc(nfft,nspden) 1693 real(dp), intent(inout) :: vtrial(nfft,nspden) 1694 1695 !Local variables------------------------------- 1696 !scalars 1697 integer :: ifft,isp 1698 real(dp) :: eigen_scf,factor,fix_resid,resid_new,resid_old 1699 character(len=500) :: message 1700 1701 ! ************************************************************************* 1702 1703 errid = AB7_NO_ERROR 1704 1705 if(nspden==4)then 1706 errid = AB7_ERROR_MIXING_ARG 1707 write(errmess, *) ' scfeig: does not work yet for nspden=4' 1708 return 1709 end if 1710 1711 !Set a fixed residual square for normalization of eigenvectors 1712 fix_resid=1.0d-4 1713 1714 !A few initialisations for the first istep 1715 if(istep==1)then 1716 1717 write(message, '(a,es12.4,a,a,a,a,a,a,a)' )& 1718 & ' scfeig: fixed PC_residual square =',fix_resid,ch10,& 1719 & ' Note that fixed resid should always be much larger',ch10,& 1720 & ' than initial PC resid square, still sufficiently',ch10,& 1721 & ' small to reduce anharmonic effects ',ch10 1722 call wrtout(std_out,message,'COLL') 1723 1724 ! Compute the preconditioned residual 1725 resid_old=0.0_dp 1726 do isp=1,nspden 1727 do ifft=1,nfft 1728 resid_old=resid_old+vrespc(ifft,isp)**2 1729 end do 1730 end do 1731 write(message, '(a,es12.4)' )' scfeig: initial PC_residual square =',resid_old 1732 call wrtout(std_out,message,'COLL') 1733 if(resid_old>1.0d-8)then 1734 errid = AB7_ERROR_MIXING_ARG 1735 write(errmess,'(a,a,a,a,a,a,a,a,a,a)') ch10,& 1736 & ' scfeig : ERROR -',ch10,& 1737 & ' This value is not good enough to allow',ch10,& 1738 & ' the computation of the eigenvectors of the SCF cycle.',ch10,& 1739 & ' It should be better than 1.0d-8 .',ch10,& 1740 & ' Action : improve the accuracy of your starting wavefunctions.' 1741 return 1742 end if 1743 1744 ! Also transfer vtrial in vtrial_old 1745 vtrial0(:,:)=vtrial(:,:) 1746 1747 ! In order to start the search for eigenvectors, 1748 ! use the tiny residual vector, renormalized 1749 factor=sqrt(fix_resid/resid_old) 1750 work(:,:,1)=vrespc(:,:)*factor 1751 vtrial(:,:)=vtrial0(:,:)+work(:,:,1) 1752 1753 ! If istep is not equal to 1 1754 else if(istep>=2)then 1755 ! 1756 ! Compute the corresponding operator expectation value 1757 ! And put the residual vector minus the difference 1758 ! between vtrial and vtrial_old 1759 ! (this is actually the action of the operator !) in vect(*,2) 1760 eigen_scf=0.0_dp 1761 do isp=1,nspden 1762 do ifft=1,nfft 1763 eigen_scf=eigen_scf+& 1764 & work(ifft,isp,1) * vrespc(ifft,isp) 1765 end do 1766 end do 1767 1768 do isp=1,nspden 1769 do ifft=1,nfft 1770 vrespc(ifft,isp)=vrespc(ifft,isp)& 1771 & +vtrial(ifft,isp)-vtrial0(ifft,isp) 1772 work(ifft,isp,2)=vrespc(ifft,isp) 1773 end do 1774 end do 1775 eigen_scf=eigen_scf/fix_resid 1776 write(message, '(a,es12.4,a)' ) & 1777 & ' scfeig : Operator expectation value ',eigen_scf,' (extremal eigenvalue * diemix)' 1778 call wrtout(std_out,message,'COLL') 1779 call wrtout(ab_out,message,'COLL') 1780 ! 1781 ! Compute residual of vect(*,2) 1782 resid_new=zero 1783 do isp=1,min(nspden,2) 1784 do ifft=1,nfft 1785 resid_new=resid_new+ work(ifft,isp,2) ** 2 1786 end do 1787 end do 1788 if (nspden==4) then 1789 do ifft=1,nfft 1790 resid_new=resid_new+two*(work(ifft,3,2)**2+work(ifft,4,2)**2) 1791 end do 1792 end if 1793 factor=sqrt(fix_resid/resid_new) 1794 if(eigen_scf<zero) then 1795 factor=-factor ! the new vector MAY be oposite to the old one 1796 ! if(factor<-one) factor=-factor ! the new vector is not opposed to the old one 1797 end if 1798 write(message, '(a,es12.4)' ) & 1799 & ' scfeig : Inverse of renormalization factor ',one/factor 1800 call wrtout(std_out,message,'COLL') 1801 call wrtout(ab_out,message,'COLL') 1802 write(message, '(a,es12.4)' ) & 1803 & ' scfeig : Convergence criterion value (->0 at convergency) ',one/factor-eigen_scf-one 1804 call wrtout(std_out,message,'COLL') 1805 call wrtout(ab_out,message,'COLL') 1806 1807 work(:,:,1)=work(:,:,2)*factor 1808 vtrial(:,:)=vtrial0(:,:)+work(:,:,1) 1809 ! End the different istep cases 1810 end if 1811 1812 end subroutine scfeig
ABINIT/sqnormm_v [ Functions ]
NAME
sqnormm_v
FUNCTION
For a series of potentials, compute square of the norm (integral over FFT grid), to obtain a square residual-like quantity (so the sum of product of values is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume). Take into account the spin components of the density and potentials (nspden), and sum over them. Need the index of the first potential to be treated, in the provided array of potentials, and the number of potentials to be treated. Might be used to compute just one square of norm, in a big array, such as to avoid copying a potential from a big array to a temporary place.
INPUTS
cplex=if 1, real space function on FFT grid is REAL, if 2, COMPLEX index=index of the first potential to be treated mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT mult=number of potentials to be treated nfft= (effective) number of FFT grid points (for this processor) npot= third dimension of the potarr array nspden=number of spin-density components opt_storage: 0, if potential is stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn] 1, if potential is stored as V, B_x, B_y, Bz (B=magn. field) potarr(cplex*nfft,nspden,npot)=array of real space potentials on FFT grid
OUTPUT
norm2(mult)= value of the square of the norm of the different potentials
SOURCE
2833 subroutine sqnormm_v(cplex,index,mpicomm, mpi_summarize,mult,nfft,norm2,npot,nspden,opt_storage,potarr) 2834 2835 !Arguments ------------------------------------ 2836 !scalars 2837 integer,intent(in) :: cplex,index,mult,nfft,npot,nspden,opt_storage,mpicomm 2838 logical, intent(in) :: mpi_summarize 2839 !arrays 2840 real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot) 2841 real(dp),intent(out) :: norm2(mult) 2842 2843 !Local variables------------------------------- 2844 !scalars 2845 integer :: ierr,ifft,ii,ispden 2846 real(dp) :: ar 2847 !arrays 2848 real(dp) :: tsec(2) 2849 2850 ! ************************************************************************* 2851 2852 !Real or complex inputs are coded 2853 DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex") 2854 DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden") 2855 2856 DBG_CHECK(index>=1,"wrong index") 2857 DBG_CHECK(mult>=1,"wrong mult") 2858 DBG_CHECK(npot>=1,"wrong npot") 2859 2860 DBG_CHECK(npot-index-mult>=-1,'npot-index-mult') 2861 2862 do ii=1,mult 2863 ar=zero 2864 do ispden=1,min(nspden,2) 2865 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar) 2866 do ifft=1,cplex*nfft 2867 ar=ar + potarr(ifft,ispden,index+ii-1)**2 2868 end do 2869 end do 2870 norm2(ii)=ar 2871 if (nspden==4) then 2872 ar=zero 2873 do ispden=3,4 2874 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar) 2875 do ifft=1,cplex*nfft 2876 ar=ar + potarr(ifft,ispden,index+ii-1)**2 2877 end do 2878 end do 2879 if (opt_storage==0) then 2880 if (cplex==1) then 2881 norm2(ii)=norm2(ii)+two*ar 2882 else 2883 norm2(ii)=norm2(ii)+ar 2884 end if 2885 else 2886 norm2(ii)=half*(norm2(ii)+ar) 2887 end if 2888 end if 2889 end do 2890 2891 !XG030513 : MPIWF reduction (addition) on norm2 is needed here 2892 if (mpi_summarize) then 2893 call timab(48,1,tsec) 2894 call xmpi_sum(norm2,mpicomm ,ierr) 2895 call timab(48,2,tsec) 2896 end if 2897 2898 end subroutine sqnormm_v
m_ab7_mixing/ab7_mixing_copy_current_step [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_copy_current_step
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
415 subroutine ab7_mixing_copy_current_step(mix, arr_resid, errid, errmess, & 416 & arr_respc, arr_paw_resid, arr_paw_respc, arr_atm) 417 418 !Arguments ------------------------------------ 419 !scalars 420 type(ab7_mixing_object), intent(inout) :: mix 421 real(dp), intent(in) :: arr_resid(mix%space * mix%nfft, mix%nspden) 422 integer, intent(out) :: errid 423 character(len = 500), intent(out) :: errmess 424 real(dp), intent(in), optional :: arr_respc(mix%space * mix%nfft, mix%nspden) 425 real(dp), intent(in), optional :: arr_paw_resid(mix%n_pawmix), arr_paw_respc(mix%n_pawmix) 426 real(dp), intent(in), optional :: arr_atm(3, mix%n_atom) 427 ! ************************************************************************* 428 429 430 if (mix%n_fftgr>0 .and. (.not. associated(mix%f_fftgr))) then 431 errid = AB7_ERROR_MIXING_ARG 432 write(errmess, '(a,a,a,a)' )ch10,& 433 & ' ab7_mixing_set_arr_current_step: ERROR (1) -',ch10,& 434 & ' Working arrays not yet allocated.' 435 return 436 end if 437 if (mix%n_pawmix>0 .and. (.not. associated(mix%f_paw))) then 438 errid = AB7_ERROR_MIXING_ARG 439 write(errmess, '(a,a,a,a)' )ch10,& 440 & ' ab7_mixing_set_arr_current_step: ERROR (2) -',ch10,& 441 & ' Working arrays not yet allocated.' 442 return 443 end if 444 if (mix%n_atom>0 .and. (.not. associated(mix%f_atm))) then 445 errid = AB7_ERROR_MIXING_ARG 446 write(errmess, '(a,a,a,a)' )ch10,& 447 & ' ab7_mixing_set_arr_current_step: ERROR (3) -',ch10,& 448 & ' Working arrays not yet allocated.' 449 return 450 end if 451 errid = AB7_NO_ERROR 452 453 if (mix%n_fftgr>0) then 454 if (mix%i_vresid(1)>0) mix%f_fftgr(:,:,mix%i_vresid(1)) = arr_resid(:,:) 455 if (present(arr_respc).and.mix%i_vrespc(1)>0) mix%f_fftgr(:,:,mix%i_vrespc(1)) = arr_respc(:,:) 456 end if 457 if (mix%n_pawmix>0) then 458 if (present(arr_paw_resid).and.mix%i_vresid(1)>0) mix%f_paw(:, mix%i_vresid(1)) = arr_paw_resid(:) 459 if (present(arr_paw_respc).and.mix%i_vrespc(1)>0) mix%f_paw(:, mix%i_vrespc(1)) = arr_paw_respc(:) 460 end if 461 if (mix%n_atom>0) then 462 if (present(arr_atm).and.mix%i_vresid(1)>0) mix%f_atm(:,:, mix%i_vresid(1)) = arr_atm(:,:) 463 end if 464 465 end subroutine ab7_mixing_copy_current_step
m_ab7_mixing/ab7_mixing_deallocate [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_deallocate
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
780 subroutine ab7_mixing_deallocate(mix) 781 782 !Arguments ------------------------------------ 783 !scalars 784 type(ab7_mixing_object), intent(inout) :: mix 785 786 !Local variables------------------------------- 787 !scalars 788 character(len = *), parameter :: subname = "ab7_mixing_deallocate" 789 ! ************************************************************************* 790 791 ABI_SFREE_PTR(mix%i_rhor) 792 ABI_SFREE_PTR(mix%i_vtrial) 793 ABI_SFREE_PTR(mix%i_vresid) 794 ABI_SFREE_PTR(mix%i_vrespc) 795 ABI_SFREE_PTR(mix%f_fftgr) 796 ABI_SFREE_PTR(mix%f_paw) 797 ABI_SFREE_PTR(mix%f_atm) 798 799 call nullify_(mix) 800 801 end subroutine ab7_mixing_deallocate
m_ab7_mixing/ab7_mixing_eval [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_eval
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
623 subroutine ab7_mixing_eval(mix, arr, istep, nfftot, ucvol, & 624 & mpi_comm, mpi_summarize, errid, errmess, & 625 & reset, isecur, pawarr, pawopt, response, etotal, potden, & 626 & resnrm, comm_atom) 627 628 !Arguments ------------------------------------ 629 !scalars 630 type(ab7_mixing_object), intent(inout) :: mix 631 integer, intent(in) :: istep, nfftot, mpi_comm 632 real(dp), intent(in) :: ucvol 633 real(dp), intent(inout) :: arr(mix%space * mix%nfft,mix%nspden) 634 logical, intent(in) :: mpi_summarize 635 integer, intent(out) :: errid 636 character(len = 500), intent(out) :: errmess 637 logical, intent(in), optional :: reset 638 integer, intent(in), optional :: isecur, comm_atom, pawopt, response 639 real(dp), intent(inout), optional, target :: pawarr(mix%n_pawmix) 640 real(dp), intent(in), optional :: etotal 641 real(dp), intent(in), optional :: potden(mix%space * mix%nfft,mix%nspden) 642 real(dp), intent(out), optional :: resnrm 643 644 !Local variables------------------------------- 645 !scalars 646 integer :: moveAtm, dbl_nnsclo, initialized, isecur_ 647 integer :: usepaw, pawoptmix_, response_ 648 real(dp) :: resnrm_ 649 !arrays 650 real(dp),target :: dum(1) 651 real(dp),pointer :: pawarr_(:) 652 653 ! ************************************************************************* 654 655 ! Argument checkings. 656 !if (mix%iscf == AB7_MIXING_NONE) then 657 ! errid = AB7_ERROR_MIXING_ARG 658 ! write(errmess, '(a,a,a,a)' )ch10,& 659 ! & ' ab7_mixing_eval: ERROR -',ch10,& 660 ! & ' No method has been chosen.' 661 ! return 662 !end if 663 if (mix%n_pawmix > 0 .and. .not. present(pawarr)) then 664 errid = AB7_ERROR_MIXING_ARG 665 write(errmess, '(a,a,a,a)' )ch10,& 666 & ' ab7_mixing_eval: ERROR -',ch10,& 667 & ' PAW is used, but no pawarr argument provided.' 668 return 669 end if 670 if (mix%n_atom > 0 .and. (.not. associated(mix%dtn_pc) .or. .not. associated(mix%xred))) then 671 errid = AB7_ERROR_MIXING_ARG 672 write(errmess, '(a,a,a,a)' )ch10,& 673 & ' ab7_mixing_eval: ERROR -',ch10,& 674 & ' Moving atoms is used, but no xred or dtn_pc attributes provided.' 675 return 676 end if 677 errid = AB7_NO_ERROR 678 679 ! Reset if requested 680 initialized = 1 681 if (present(reset)) then 682 if (reset) initialized = 0 683 end if 684 685 ! Miscellaneous 686 moveAtm = 0 687 if (mix%n_atom > 0) moveAtm = 1 688 isecur_ = 0 689 if (present(isecur)) isecur_ = isecur 690 usepaw = 0 691 if (mix%n_pawmix > 0) usepaw = 1 692 pawoptmix_ = 0 693 if (present(pawopt)) pawoptmix_ = pawopt 694 response_ = 0 695 if (present(response)) response_ = response 696 pawarr_ => dum ; if (present(pawarr)) pawarr_ => pawarr 697 698 ! Do the mixing. 699 resnrm_ = 0.d0 700 if (mix%iscf == AB7_MIXING_NONE) then 701 arr(:,:)=arr(:,:)+mix%f_fftgr(:,:,1) 702 else if (mix%iscf == AB7_MIXING_EIG) then 703 ! This routine compute the eigenvalues of the SCF operator 704 call scfeig(istep, mix%space * mix%nfft, mix%nspden, & 705 & mix%f_fftgr(:,:,mix%i_vrespc(1)), arr, & 706 & mix%f_fftgr(:,:,1), mix%f_fftgr(:,:,4:5), errid, errmess) 707 else if (mix%iscf == AB7_MIXING_SIMPLE .or. & 708 & mix%iscf == AB7_MIXING_ANDERSON .or. & 709 & mix%iscf == AB7_MIXING_ANDERSON_2 .or. & 710 & mix%iscf == AB7_MIXING_PULAY) then 711 if (present(comm_atom)) then 712 call scfopt(mix%space, mix%f_fftgr,mix%f_paw,mix%iscf,istep,& 713 & mix%i_vrespc,mix%i_vtrial, & 714 & mpi_comm,mpi_summarize,mix%nfft,mix%n_pawmix,mix%nspden, & 715 & mix%n_fftgr,mix%n_index,mix%kind,pawoptmix_,usepaw,pawarr_, & 716 & resnrm_, arr, errid, errmess, comm_atom=comm_atom) 717 else 718 call scfopt(mix%space, mix%f_fftgr,mix%f_paw,mix%iscf,istep,& 719 & mix%i_vrespc,mix%i_vtrial, & 720 & mpi_comm,mpi_summarize,mix%nfft,mix%n_pawmix,mix%nspden, & 721 & mix%n_fftgr,mix%n_index,mix%kind,pawoptmix_,usepaw,pawarr_, & 722 & resnrm_, arr, errid, errmess) 723 end if 724 ! Change atomic positions 725 if((istep==1 .or. mix%iscf==AB7_MIXING_SIMPLE) .and. mix%n_atom > 0)then 726 ! GAF: 2009-06-03 727 ! Apparently there are not reason 728 ! to restrict iscf=2 for ionmov=5 729 mix%xred(:,:) = mix%xred(:,:) + mix%dtn_pc(:,:) 730 end if 731 else if (mix%iscf == AB7_MIXING_CG_ENERGY .or. mix%iscf == AB7_MIXING_CG_ENERGY_2) then 732 ! Optimize next vtrial using an algorithm based 733 ! on the conjugate gradient minimization of etotal 734 if (.not. present(etotal) .or. .not. present(potden)) then 735 errid = AB7_ERROR_MIXING_ARG 736 write(errmess, '(a,a,a,a)' )ch10,& 737 & ' ab7_mixing_eval: ERROR -',ch10,& 738 & ' Arguments etotal or potden are missing for CG on energy methods.' 739 return 740 end if 741 if (mix%n_atom == 0) then 742 ABI_MALLOC(mix%xred,(3,0)) 743 ABI_MALLOC(mix%dtn_pc,(3,0)) 744 end if 745 call scfcge(mix%space,dbl_nnsclo,mix%dtn_pc,etotal,mix%f_atm,& 746 & mix%f_fftgr,initialized,mix%iscf,isecur_,istep,& 747 & mix%i_rhor,mix%i_vresid,mix%i_vrespc,moveAtm,& 748 & mpi_comm,mpi_summarize,mix%n_atom,mix%nfft,nfftot,& 749 & mix%nspden,mix%n_fftgr,mix%n_index,mix%kind,& 750 & response_,potden,ucvol,arr,mix%xred, errid, errmess) 751 if (mix%n_atom == 0) then 752 ABI_FREE(mix%xred) 753 ABI_FREE(mix%dtn_pc) 754 end if 755 if (dbl_nnsclo == 1) errid = AB7_ERROR_MIXING_INC_NNSLOOP 756 end if 757 758 if (present(resnrm)) resnrm = resnrm_ 759 760 end subroutine ab7_mixing_eval
m_ab7_mixing/ab7_mixing_eval_allocate [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_eval_allocate
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
485 subroutine ab7_mixing_eval_allocate(mix, istep) 486 487 !Arguments ------------------------------------ 488 !scalars 489 type(ab7_mixing_object), intent(inout) :: mix 490 integer, intent(in), optional :: istep 491 492 !Local variables------------------------------- 493 !scalars 494 integer :: istep_,temp_unit !, i_stat 495 real(dp) :: tsec(2) 496 character(len = *), parameter :: subname = "ab7_mixing_eval_allocate" 497 character(len=500) :: msg 498 499 ! ************************************************************************* 500 501 istep_ = 1 502 if (present(istep)) istep_ = istep 503 504 ! Allocate work array. 505 if (.not. associated(mix%f_fftgr)) then 506 !allocate(mix%f_fftgr(mix%space * mix%nfft,mix%nspden,mix%n_fftgr), stat = i_stat) 507 !call memocc_abi(i_stat, mix%f_fftgr, 'mix%f_fftgr', subname) 508 ABI_MALLOC(mix%f_fftgr,(mix%space * mix%nfft,mix%nspden,mix%n_fftgr)) 509 mix%f_fftgr(:,:,:)=zero 510 if (mix%mffmem == 0 .and. istep_ > 1 .and. mix%n_fftgr>0) then 511 call timab(83,1,tsec) 512 if (open_file(mix%diskCache,msg,newunit=temp_unit,form='unformatted',status='old') /= 0) then 513 ABI_ERROR(msg) 514 end if 515 rewind(temp_unit) 516 read(temp_unit) mix%f_fftgr 517 if (mix%n_pawmix == 0) close(unit=temp_unit) 518 call timab(83,2,tsec) 519 end if 520 end if 521 ! Allocate PAW work array. 522 if (.not. associated(mix%f_paw)) then 523 !allocate(mix%f_paw(mix%n_pawmix,mix%n_fftgr), stat = i_stat) 524 !call memocc_abi(i_stat, mix%f_paw, 'mix%f_paw', subname) 525 ABI_MALLOC(mix%f_paw,(mix%n_pawmix,mix%n_fftgr)) 526 if (mix%n_pawmix > 0 .and. mix%n_fftgr>0) then 527 mix%f_paw(:,:)=zero 528 if (mix%mffmem == 0 .and. istep_ > 1) then 529 read(temp_unit) mix%f_paw 530 close(unit=temp_unit) 531 call timab(83,2,tsec) 532 end if 533 end if 534 end if 535 ! Allocate atom work array. 536 if (.not. associated(mix%f_atm)) then 537 !allocate(mix%f_atm(3,mix%n_atom,mix%n_fftgr), stat = i_stat) 538 !call memocc_abi(i_stat, mix%f_atm, 'mix%f_atm', subname) 539 ABI_MALLOC(mix%f_atm,(3,mix%n_atom,mix%n_fftgr)) 540 end if 541 542 end subroutine ab7_mixing_eval_allocate
m_ab7_mixing/ab7_mixing_eval_deallocate [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_eval_deallocate
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
562 subroutine ab7_mixing_eval_deallocate(mix) 563 564 !Arguments ------------------------------------ 565 !scalars 566 type(ab7_mixing_object), intent(inout) :: mix 567 568 !Local variables------------------------------- 569 !scalars 570 integer :: temp_unit !i_all, i_stat 571 real(dp) :: tsec(2) 572 character(len = *), parameter :: subname = "ab7_mixing_eval_deallocate" 573 character(len=500) :: msg 574 575 ! ************************************************************************* 576 577 ! Save on disk and deallocate work array in case on disk cache only. 578 if (mix%mffmem == 0) then 579 call timab(83,1,tsec) 580 if (open_file(mix%diskCache,msg,newunit=temp_unit,form='unformatted',status='unknown') /= 0) then 581 ABI_ERROR(msg) 582 end if 583 rewind(temp_unit) 584 ! VALGRIND complains not all of f_fftgr_disk is initialized 585 if (mix%n_fftgr > 0) then 586 write(temp_unit) mix%f_fftgr 587 end if 588 if (mix%n_pawmix > 0 .and. mix%n_fftgr > 0) then 589 write(temp_unit) mix%f_paw 590 end if 591 close(unit=temp_unit) 592 call timab(83,2,tsec) 593 if (associated(mix%f_fftgr)) then 594 ABI_FREE(mix%f_fftgr) 595 nullify(mix%f_fftgr) 596 end if 597 if (associated(mix%f_paw)) then 598 ABI_FREE(mix%f_paw) 599 nullify(mix%f_paw) 600 end if 601 end if 602 603 end subroutine ab7_mixing_eval_deallocate
m_ab7_mixing/ab7_mixing_new [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_new
FUNCTION
INPUTS
OUTPUT
NOTES
SOURCE
161 subroutine ab7_mixing_new(mix, iscf, kind, space, nfft, nspden, & 162 & npawmix, errid, errmess, npulayit, useprec) 163 164 !Arguments ------------------------------------ 165 !scalars 166 type(ab7_mixing_object), intent(out) :: mix 167 integer, intent(in) :: iscf, kind, space, nfft, nspden, npawmix 168 integer, intent(out) :: errid 169 character(len = 500), intent(out) :: errmess 170 integer, intent(in), optional :: npulayit 171 logical, intent(in), optional :: useprec 172 173 !Local variables------------------------------- 174 !scalars 175 integer :: ii !, i_stat 176 character(len = *), parameter :: subname = "ab7_mixing_new" 177 ! ************************************************************************* 178 179 ! Set default values. 180 call init_(mix) 181 182 ! Argument checkings. 183 if (kind /= AB7_MIXING_POTENTIAL .and. kind /= AB7_MIXING_DENSITY) then 184 errid = AB7_ERROR_MIXING_ARG 185 write(errmess, '(a,a,a,a)' )ch10,& 186 & ' ab7_mixing_set_arrays: ERROR -',ch10,& 187 & ' Mixing must be done on density or potential only.' 188 return 189 end if 190 if (space /= AB7_MIXING_REAL_SPACE .and. space /= AB7_MIXING_FOURRIER_SPACE) then 191 errid = AB7_ERROR_MIXING_ARG 192 write(errmess, '(a,a,a,a)' )ch10,& 193 & ' ab7_mixing_set_arrays: ERROR -',ch10,& 194 & ' Mixing must be done in real or Fourrier space only.' 195 return 196 end if 197 if (iscf /= AB7_MIXING_EIG .and. iscf /= AB7_MIXING_SIMPLE .and. & 198 & iscf /= AB7_MIXING_ANDERSON .and. & 199 & iscf /= AB7_MIXING_ANDERSON_2 .and. & 200 & iscf /= AB7_MIXING_CG_ENERGY .and. & 201 & iscf /= AB7_MIXING_PULAY .and. & 202 & iscf /= AB7_MIXING_CG_ENERGY_2 .and. & 203 & iscf /= AB7_MIXING_NONE) then 204 errid = AB7_ERROR_MIXING_ARG 205 write(errmess, "(A,I0,A)") "Unknown mixing scheme (", iscf, ")." 206 return 207 end if 208 errid = AB7_NO_ERROR 209 210 ! Mandatory arguments. 211 mix%iscf = iscf 212 mix%kind = kind 213 mix%space = space 214 mix%nfft = nfft 215 mix%nspden = nspden 216 mix%n_pawmix = npawmix 217 218 ! Optional arguments. 219 if (present(useprec)) mix%useprec = useprec 220 221 ! Set-up internal dimensions. 222 !These arrays are needed only in the self-consistent case 223 if (iscf == AB7_MIXING_NONE) then 224 ! For iscf==0, one additional vector is needed. 225 ! The index 1 is attributed to the new residual potential. 226 mix%n_fftgr=1 ; mix%n_index=1 227 else if (iscf == AB7_MIXING_EIG) then 228 ! For iscf==1, five additional vectors are needed. 229 ! The index 1 is attributed to the old trial potential, 230 ! The new residual potential, and the new 231 ! preconditioned residual potential receive now a temporary index 232 ! The indices number 4 and 5 are attributed to work vectors. 233 mix%n_fftgr=5 ; mix%n_index=1 234 else if(iscf == AB7_MIXING_SIMPLE) then 235 ! For iscf==2, three additional vectors are needed. 236 ! The index number 1 is attributed to the old trial vector 237 ! The new residual potential, and the new preconditioned 238 ! residual potential, receive now a temporary index. 239 mix%n_fftgr=3 ; mix%n_index=1 240 if (.not. mix%useprec) mix%n_fftgr = 2 241 else if(iscf == AB7_MIXING_ANDERSON) then 242 ! For iscf==3 , four additional vectors are needed. 243 ! The index number 1 is attributed to the old trial vector 244 ! The new residual potential, and the new and old preconditioned 245 ! residual potential, receive now a temporary index. 246 mix%n_fftgr=4 ; mix%n_index=2 247 if (.not. mix%useprec) mix%n_fftgr = 3 248 else if (iscf == AB7_MIXING_ANDERSON_2) then 249 ! For iscf==4 , six additional vectors are needed. 250 ! The indices number 1 and 2 are attributed to two old trial vectors 251 ! The new residual potential, and the new and two old preconditioned 252 ! residual potentials, receive now a temporary index. 253 mix%n_fftgr=6 ; mix%n_index=3 254 if (.not. mix%useprec) mix%n_fftgr = 5 255 else if(iscf == AB7_MIXING_CG_ENERGY .or. iscf == AB7_MIXING_CG_ENERGY_2) then 256 ! For iscf==5 or 6, ten additional vectors are needed 257 ! The index number 1 is attributed to the old trial vector 258 ! The index number 6 is attributed to the search vector 259 ! Other indices are attributed now. Altogether ten vectors 260 mix%n_fftgr=10 ; mix%n_index=3 261 else if(iscf == AB7_MIXING_PULAY) then 262 ! For iscf==7, lot of additional vectors are needed 263 ! The index number 1 is attributed to the old trial vector 264 ! The index number 2 is attributed to the old residual 265 ! The indices number 2 and 3 are attributed to two old precond. residuals 266 ! Other indices are attributed now. 267 if (present(npulayit)) mix%n_pulayit = npulayit 268 mix%n_fftgr=2+2*mix%n_pulayit ; mix%n_index=1+mix%n_pulayit 269 if (.not. mix%useprec) mix%n_fftgr = 1+2*mix%n_pulayit 270 end if ! iscf cases 271 272 ! Allocate new arrays. 273 !allocate(mix%i_rhor(mix%n_index), stat = i_stat) 274 !call memocc_abi(i_stat, mix%i_rhor, 'mix%i_rhor', subname) 275 ABI_MALLOC(mix%i_rhor,(mix%n_index)) 276 mix%i_rhor(:)=0 277 !allocate(mix%i_vtrial(mix%n_index), stat = i_stat) 278 !call memocc_abi(i_stat, mix%i_vtrial, 'mix%i_vtrial', subname) 279 ABI_MALLOC(mix%i_vtrial,(mix%n_index)) 280 mix%i_vtrial(:)=0 281 !allocate(mix%i_vresid(mix%n_index), stat = i_stat) 282 !call memocc_abi(i_stat, mix%i_vresid, 'mix%i_vresid', subname) 283 ABI_MALLOC(mix%i_vresid,(mix%n_index)) 284 mix%i_vresid(:)=0 285 !allocate(mix%i_vrespc(mix%n_index), stat = i_stat) 286 !call memocc_abi(i_stat, mix%i_vrespc, 'mix%i_vrespc', subname) 287 ABI_MALLOC(mix%i_vrespc,(mix%n_index)) 288 mix%i_vrespc(:)=0 289 290 ! Setup initial values. 291 if (iscf == AB7_MIXING_NONE) then 292 mix%i_vresid(1)=1 293 else if (iscf == AB7_MIXING_EIG) then 294 mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 ; mix%i_vrespc(1)=3 295 else if(iscf == AB7_MIXING_SIMPLE) then 296 mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 ; mix%i_vrespc(1)=3 297 if (.not. mix%useprec) mix%i_vrespc(1)=2 298 else if(iscf == AB7_MIXING_ANDERSON) then 299 mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 300 if (mix%useprec) then 301 mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=4 302 else 303 mix%i_vrespc(1)=2 ; mix%i_vrespc(2)=3 304 end if 305 else if (iscf == AB7_MIXING_ANDERSON_2) then 306 mix%i_vtrial(1)=1 ; mix%i_vtrial(2)=2 307 mix%i_vresid(1)=3 308 if (mix%useprec) then 309 mix%i_vrespc(1)=4 ; mix%i_vrespc(2)=5 ; mix%i_vrespc(3)=6 310 else 311 mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=4 ; mix%i_vrespc(3)=5 312 end if 313 else if(iscf == AB7_MIXING_CG_ENERGY .or. iscf == AB7_MIXING_CG_ENERGY_2) then 314 mix%n_fftgr=10 ; mix%n_index=3 315 mix%i_vtrial(1)=1 316 mix%i_vresid(1)=2 ; mix%i_vresid(2)=4 ; mix%i_vresid(3)=7 317 mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=5 ; mix%i_vrespc(3)=8 318 mix%i_rhor(2)=9 ; mix%i_rhor(3)=10 319 else if(iscf == AB7_MIXING_PULAY) then 320 do ii=1,mix%n_pulayit 321 mix%i_vtrial(ii)=2*ii-1 ; mix%i_vrespc(ii)=2*ii 322 end do 323 mix%i_vrespc(mix%n_pulayit+1)=2*mix%n_pulayit+1 324 mix%i_vresid(1)=2*mix%n_pulayit+2 325 if (.not. mix%useprec) mix%i_vresid(1)=2 326 end if ! iscf cases 327 328 end subroutine ab7_mixing_new
m_ab7_mixing/ab7_mixing_use_disk_cache [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_use_disk_cache
FUNCTION
INPUTS
OUTPUT
NOTES
Obsolete?
SOURCE
346 subroutine ab7_mixing_use_disk_cache(mix, fnametmp_fft) 347 348 !Arguments ------------------------------------ 349 !scalars 350 type(ab7_mixing_object), intent(inout) :: mix 351 character(len = *), intent(in) :: fnametmp_fft 352 ! ************************************************************************* 353 354 if (len(trim(fnametmp_fft)) > 0) then 355 mix%mffmem = 0 356 write(mix%diskCache, "(A)") fnametmp_fft 357 else 358 mix%mffmem = 1 359 end if 360 361 end subroutine ab7_mixing_use_disk_cache
m_ab7_mixing/ab7_mixing_use_moving_atoms [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
ab7_mixing_use_moving_atoms
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
381 subroutine ab7_mixing_use_moving_atoms(mix, natom, xred, dtn_pc) 382 383 !Arguments ------------------------------------ 384 !scalars 385 type(ab7_mixing_object), intent(inout) :: mix 386 integer, intent(in) :: natom 387 real(dp), intent(in), target :: dtn_pc(3, natom) 388 real(dp), intent(in), target :: xred(3, natom) 389 390 ! ************************************************************************* 391 392 mix%n_atom = natom 393 mix%dtn_pc => dtn_pc 394 mix%xred => xred 395 396 end subroutine ab7_mixing_use_moving_atoms
m_ab7_mixing/init_ [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
init_
FUNCTION
Initialize the object
SOURCE
96 subroutine init_(mix) 97 98 !Arguments ------------------------------------ 99 !scalars 100 type(ab7_mixing_object), intent(out) :: mix 101 ! ************************************************************************* 102 103 ! Default values. 104 mix%iscf = AB7_MIXING_NONE 105 mix%mffmem = 1 106 mix%n_index = 0 107 mix%n_fftgr = 0 108 mix%n_pulayit = 7 109 mix%n_pawmix = 0 110 mix%n_atom = 0 111 mix%space = 0 112 mix%useprec = .true. 113 114 call nullify_(mix) 115 116 end subroutine init_
m_ab7_mixing/nullify [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
nullify_
FUNCTION
Nullify the pointers
SOURCE
128 subroutine nullify_(mix) 129 130 !Arguments ------------------------------------ 131 !scalars 132 type(ab7_mixing_object), intent(inout) :: mix 133 ! ************************************************************************* 134 135 ! Nullify internal pointers. 136 nullify(mix%i_rhor) 137 nullify(mix%i_vtrial) 138 nullify(mix%i_vresid) 139 nullify(mix%i_vrespc) 140 nullify(mix%f_fftgr) 141 nullify(mix%f_atm) 142 nullify(mix%f_paw) 143 144 end subroutine nullify_
m_ab7_mixing/scfcge [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
scfcge
FUNCTION
Compute the next vtrial of the SCF cycle. Uses a conjugate gradient minimization of the total energy Can move only the trial potential (if moved_atm_inside==0), or move the trial atomic positions as well (if moved_atm_inside==1).
INPUTS
cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX dtn_pc(3,natom)=preconditioned change of atomic position, in reduced coordinates. Will be quickly transferred to f_atm(:,:,i_vrespc(1)) etotal=the actual total energy initialized= if 0, the initialization of the gstate run is not yet finished iscf =5 => SCF cycle, CG based on estimation of energy gradient =6 => SCF cycle, CG based on true minimization of the energy isecur=level of security of the computation istep= number of the step in the SCF cycle moved_atm_inside: if==1, the atoms are allowed to move. mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT natom=number of atoms nfft=(effective) number of FFT grid points (for this processor) nfftot=total number of FFT grid points nspden=number of spin-density components n_fftgr=third dimension of the array f_fftgr n_index=dimension for indices of potential/density (see i_vresid, ivrespc, i_rhor...) opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential 1 vtrial (and also f_fftgr) actually contains the trial density response= if 0, GS calculation, if 1, RF calculation, intrinsically harmonic ! rhor(cplex*nfft,nspden)=actual density ucvol=unit cell volume in bohr**3
OUTPUT
dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
SIDE EFFECTS
Input/Output:
vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
the input residual of the potential and Hellman-Feynman forces
at output, it is the new trial potential .
xred(3,natom)=(needed if moved_atm_inside==1)
reduced dimensionless atomic coordinates
at input, those that generated the input residual of the potential
and Hellman-Feynman forces, at output, these are the new ones.
f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
The input vtrial is transferred, at output, in f_fftgr(:,:,1).
The input f_fftgr(:,:,i_vresid(1)) contains the last residual.
the value of i_vresid(1) is transferred to i_vresid(2) at output.
The input f_fftgr(:,:,i_vresid(2)) contains the old residual.
the value of i_vresid(2) is transferred to i_vresid(3) at output.
The input f_fftgr(:,:,i_vresid(3)) contains the previous last residual.
For the preconditioned potential residual, the same logic as for the
the potential residual is used, with i_vrespc replacing i_vresid.
The input rhor is transferred, at output, in f_fft(:,:,i_rhor(2)).
The old density is input in f_fft(:,:,i_rhor(2)), and the value of
i_rhor(2) is transferred to i_rhor(3) before the end of the routine.
The input/output search vector is stored in f_fftgr(:,:,6)
f_atm(3,natom,n_fftgr)=different functions defined for each atom :
The input xred is transferred, at output, in f_atm(:,:,1).
The input f_atm(:,:,i_vresid(1)) contains minus the HF forces.
the value of i_vresid(1) is transferred to i_vresid(2) at output.
The input f_atm(:,:,i_vresid(2)) contains minus the old HF forces.
the value of i_vresid(2) is transferred to i_vresid(3) at output.
The input f_atm(:,:,i_vresid(3)) contains minus the previous old HF forces.
For the preconditioned change of atomic positions, the same logic as for the
the potential residual is used, with i_vrespc replacing i_vresid.
The input/output search vector is stored in f_atm(:,:,6)
i_rhor(2:3)=index of the density (past and previous past) in the array f_fftgr
i_vresid(3)=index of the residual potentials (present, past and previous
past) in the array f_fftgr; also similar index for minus Hellman-Feynman
forces in the array f_atm .
i_vrespc(3)=index of the preconditioned residual potentials
(present, past and previous past) in the array f_fftgr ;
also similar index for the preconditioned change of atomic position (dtn_pc).
TODO
This routine is much too difficult to read ! Should be rewritten ... Maybe make separate subroutines for line search and CG step ?!
SOURCE
889 subroutine scfcge(cplex,dbl_nnsclo,dtn_pc,etotal,f_atm,& 890 & f_fftgr,initialized,iscf,isecur,istep,& 891 & i_rhor,i_vresid,i_vrespc,moved_atm_inside,mpicomm,mpi_summarize,& 892 & natom,nfft,nfftot,nspden,n_fftgr,n_index,opt_denpot,response,rhor,ucvol,vtrial,xred,errid,errmess) 893 894 !Arguments ------------------------------------ 895 !scalars 896 integer,intent(in) :: cplex,initialized,iscf,isecur,istep,moved_atm_inside,mpicomm 897 integer,intent(in) :: n_fftgr,n_index,natom,nfft,nfftot,nspden,opt_denpot,response 898 integer,intent(out) :: dbl_nnsclo, errid 899 character(len = 500), intent(out) :: errmess 900 logical, intent(in) :: mpi_summarize 901 real(dp),intent(in) :: etotal,ucvol 902 !arrays 903 integer,intent(inout) :: i_rhor(n_index),i_vresid(n_index),i_vrespc(n_index) 904 real(dp),intent(in) :: dtn_pc(3,natom),rhor(cplex*nfft,nspden) 905 real(dp),intent(inout) :: f_atm(3,natom,n_fftgr) 906 real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr) 907 real(dp),intent(inout) :: vtrial(cplex*nfft,nspden),xred(3,natom) 908 909 !Local variables------------------------------- 910 !mlinmin gives the maximum number of steps in the line minimization 911 ! after which the algorithm is restarted (with a decrease of the 912 ! adaptative trial step length). This number should not be large, 913 ! since if the potential landscape is harmonic, the number of 914 ! search steps should be small. If it is large, we are not in the 915 ! harmonic region, and the CG algorithm will not be really useful, 916 ! so one can just restart the algorithm ... 917 !scalars 918 integer,parameter :: mlinmin=5 919 integer,save :: end_linmin,iline_cge,ilinear,ilinmin,isecur_eff,nlinear 920 integer,save :: number_of_restart,status 921 integer :: choice,iatom,idir,ifft,iline_cge_input,ilinmin_input,isp 922 integer :: testcg,tmp,errid_ 923 real(dp),save :: d2edv2_old2,d_lambda_old2,dedv_old2,etotal_old 924 real(dp),save :: etotal_previous=MAGIC_UNDEF,lambda_adapt,lambda_new,lambda_old,resid_old 925 real(dp) :: d2e11,d2e12,d2e22,d2edv2_new,d2edv2_old 926 real(dp) :: d2edv2_predict,d_lambda,de1,de2,dedv_mix 927 real(dp) :: dedv_new,dedv_old,dedv_predict,determ,etotal_input 928 real(dp) :: etotal_predict,gamma,lambda_input,lambda_predict2 929 real(dp) :: lambda_predict=1.0_dp,ratio,reduction 930 real(dp) :: resid_input,temp 931 character(len=500) :: message 932 !arrays 933 real(dp) :: resid_new(1) 934 real(dp), allocatable :: tmp_fft1(:,:) 935 936 ! ************************************************************************* 937 938 errid = AB7_NO_ERROR 939 dbl_nnsclo = 0 940 941 !reduction gives the level of reduction of the error in 942 !the line minimization to be reached for the minimization to be 943 !considered successfull 944 reduction=0.1_dp 945 946 !nlinear increases with the number of times the 2D minimization succeded 947 !to reach the true minimum directly. It is a measure of the 948 !degree of parabolicity of the problem, and is used to 949 !skip some steps by performing extrapolation. 950 if(istep==1)then 951 952 ! Skipping some steps is sometimes unsecure, so it is possible 953 ! to make nlinear start at a negative value - if isecur is positive 954 isecur_eff=isecur 955 nlinear=min(-isecur_eff,0) 956 ilinear=0 957 958 ! Response function calculation are intrinsically harmonic, so one 959 ! can shift isecur (by -2), and start with a positive nlinear 960 if(response==1)then 961 isecur_eff=isecur-2 962 nlinear=-isecur_eff 963 ilinear=nlinear 964 end if 965 966 iline_cge=0 967 ilinmin=0 968 end if 969 970 !Compute actual residual resid_new (residual of f_fftgr(:,:,i_vrespc(1)) 971 call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr) 972 973 !Save input residual and ilinmin for final printing 974 resid_input=resid_new(1) 975 etotal_input=etotal 976 ilinmin_input=ilinmin 977 iline_cge_input=iline_cge 978 !Transfer dtn_pc in f_atm 979 if(moved_atm_inside==1)then 980 f_atm(:,:,i_vrespc(1))=dtn_pc(:,:) 981 end if 982 983 !======================================================================= 984 !Now the routine is decomposed in three mutually exclusive parts : 985 !if(istep==1)then initialize the algorithm 986 !else if(ilinmin>0)then perform the line minimisation 987 !else if(ilinmin==0)then determine the new search direction (CG step) 988 !======================================================================= 989 990 991 !-------------------------------------- 992 !Here initialize the algorithm 993 if(istep==1)then 994 995 ! At the beginning of each gstate run, lambda_adapt is forced to have the 996 ! same value, that is 1.0_dp. In the other cases when istep=1 (at different 997 ! broyden steps, for example), the previously obtained 998 ! adaptive value is kept. 999 if(initialized==0)lambda_adapt=1.0_dp 1000 lambda_old=0.0_dp 1001 lambda_input=0.0_dp 1002 number_of_restart=0 1003 lambda_new=lambda_adapt 1004 1005 f_fftgr(:,:,1)=vtrial(:,:) 1006 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 1007 1008 ! This copy must be written in F77, because of stack problems on the DECs 1009 do isp=1,nspden 1010 do ifft=1,cplex*nfft 1011 f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1)) 1012 end do 1013 end do 1014 vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6) 1015 if(moved_atm_inside==1)then 1016 f_atm(:,:,1)=xred(:,:) 1017 f_atm(:,:,i_rhor(2))=xred(:,:) 1018 ! There shouldn t be problems with the stack size for this small array. 1019 f_atm(:,:,6)=f_atm(:,:,i_vrespc(1)) 1020 xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6) 1021 end if 1022 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 1023 tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 1024 ilinmin=1 1025 resid_old=resid_new(1) 1026 etotal_old=etotal 1027 1028 status=0 1029 1030 ! -------------------------------------- 1031 1032 ! Here performs the line minimisation 1033 else if(ilinmin>0)then 1034 1035 lambda_input=lambda_new 1036 1037 ! The choice with the Brent algorithm has been abandoned in version 1.6.m 1038 1039 ! Compute the approximate energy derivatives dedv_new and dedv_old, 1040 ! from vresid and vresid_old 1041 choice=2 1042 call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,& 1043 & f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,& 1044 & natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred) 1045 d_lambda=lambda_new-lambda_old 1046 dedv_old=dedv_old/d_lambda 1047 dedv_new=dedv_new/d_lambda 1048 1049 ! DEBUG 1050 ! write(std_out,'(a,4es12.4,i3)' )' scfcge:lold,lnew,dold,dnew,status', & 1051 ! & lambda_old,lambda_new,dedv_old,dedv_new,status 1052 ! ENDDEBUG 1053 1054 if(status==0 .or. status==3)then 1055 ! 1056 ! Then, compute a predicted point along the line 1057 ! The value of choice determines the minimization algorithm 1058 ! choice=1 uses the two values of the derivative of the energy 1059 ! choice=2 uses the two values of the energy, and and estimate of the 1060 ! second derivative at the mid-point. 1061 1062 choice=1 1063 if(iscf==6)choice=2 1064 call findminscf(choice,dedv_new,dedv_old,dedv_predict,& 1065 & d2edv2_new,d2edv2_old,d2edv2_predict,& 1066 & etotal,etotal_old,etotal_predict,& 1067 & lambda_new,lambda_old,lambda_predict,errid_,message) 1068 if (errid_ /= AB7_NO_ERROR) then 1069 call wrtout(std_out,message,'COLL') 1070 end if 1071 1072 ! Suppress the next line for debugging (there is another such line) 1073 status=0 1074 1075 ! DEBUG 1076 ! Keep this debugging feature : it gives access to the investigation of lines 1077 ! in a different approach 1078 ! if(response==1 .and. istep>8)then 1079 ! lambda_predict=1.2d-2 1080 ! if(istep>=15)lambda_predict=lambda_predict-0.002 1081 ! if(istep>=14)stop 1082 ! status=3 1083 ! end if 1084 ! ENDDEBUG 1085 1086 else 1087 if(status/=-1)then 1088 status=-1 1089 lambda_predict=-2.5_dp 1090 else 1091 lambda_predict=lambda_predict+0.1_dp 1092 end if 1093 end if 1094 1095 ! If the predicted point is very close to the most recent 1096 ! computed point, while this is the first trial on this line, 1097 ! then we are in the linear regime : 1098 ! nlinear is increased by one unit. For the time being, do this even when 1099 ! moved_atm_inside==1 (the code still works when it is done, but it 1100 ! seems to be a bit unstable). The maximal value of nlinear is 1, except 1101 ! when isecur_eff is a negative number, less than -1. 1102 if( abs(lambda_predict-lambda_new)/& 1103 & (abs(lambda_predict)+abs(lambda_new)) < 0.01 .and. ilinmin==1 ) then 1104 ! if(moved_atm_inside==0 .and. nlinear<max(1,-isecur_eff) )nlinear=nlinear+1 1105 if(nlinear<max(1,-isecur_eff))nlinear=nlinear+1 1106 ilinear=nlinear 1107 end if 1108 1109 ! If the predicted point is close to the most recent computed point, 1110 ! or the previous one, set on the flag of end of line minization 1111 end_linmin=0 1112 if(abs(lambda_new-lambda_predict)*2.0_dp& 1113 & /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1 1114 if(abs(lambda_old-lambda_predict)*2.0_dp& 1115 & /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1 1116 1117 if(status/=0)end_linmin=0 1118 1119 ! Save the closest old lambda, if needed, 1120 ! also examine the reduction of the interval, and eventual stop 1121 ! the present line minimisation, because of convergence (end_linmin=1) 1122 ! Also treat the case in which the predicted value of lambda is negative, 1123 ! or definitely too small in which case the algorithm has to be restarted 1124 ! (not a very good solution, though ...) 1125 ! Finally also treat the case where insufficiently converged 1126 ! density at lambda=0.0_dp happens, which screws up the line minimisation. 1127 1128 ! Here restart the algorithm with the best vtrial. 1129 ! Also make reduction in lambda_adapt 1130 ! DEBUG 1131 ! write(std_out,*)' scfcge : status=',status 1132 ! ENDDEBUG 1133 if( end_linmin==0 .and. status==0 .and. & 1134 & ( (lambda_predict<0.005_dp*lambda_adapt .and. iscf==5) .or. & 1135 & (abs(lambda_predict)<0.005_dp*lambda_adapt .and. iscf==6).or. & 1136 & ilinmin==mlinmin ) )then 1137 if(number_of_restart>12)then 1138 errid = AB7_ERROR_MIXING_CONVERGENCE 1139 write(errmess,'(a,a,i0,a,a,a,a,a)')& 1140 & 'Potential-based CG line minimization not',' converged after ',number_of_restart,' restarts. ',ch10,& 1141 & 'Action : read the eventual warnings about lack of convergence.',ch10,& 1142 & 'Some might be relevant. Otherwise, raise nband. Returning' 1143 ABI_WARNING(errmess) 1144 return 1145 end if 1146 ! Make reduction in lambda_adapt (kind of steepest descent...) 1147 write(message,'(a,a,a)')& 1148 & 'Potential-based CG line minimization has trouble to converge.',ch10,& 1149 & 'The algorithm is restarted with more secure parameters.' 1150 ABI_WARNING(message) 1151 number_of_restart=number_of_restart+1 1152 ! At the second restart, double the number of non-self consistent loops. 1153 if(number_of_restart>=2)dbl_nnsclo=1 1154 lambda_adapt=lambda_adapt*0.7_dp 1155 lambda_new=lambda_adapt 1156 ! If the last energy is better than the old one, transfer the data. 1157 ! Otherwise, no transfer must occur (very simple to code...) 1158 if(etotal<etotal_old .or. abs(lambda_old)<1.0d-8)then 1159 f_fftgr(:,:,1)=vtrial(:,:) 1160 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 1161 do isp=1,nspden 1162 do ifft=1,cplex*nfft 1163 f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1)) 1164 end do 1165 end do 1166 if(moved_atm_inside==1)then 1167 f_atm(:,:,1)=xred(:,:) 1168 f_atm(:,:,i_rhor(2))=xred(:,:) 1169 f_atm(:,:,6)=f_atm(:,:,i_vrespc(1)) 1170 end if 1171 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 1172 tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 1173 resid_old=resid_new(1) 1174 etotal_old=etotal 1175 end if 1176 lambda_old=0.0_dp 1177 ilinmin=1 1178 ! Putting the flag to -1 avoids the usual actions taken with end_linmin=1 1179 end_linmin=-1 1180 ! Also put ilinear and nlinear to 0 1181 ilinear=0 1182 nlinear=0 1183 1184 ! Here lambda_new is the closest to lambda_predict, 1185 ! or lambda_old is still 0.0_dp, while the energy shows that the minimum 1186 ! is away from 0.0_dp (insufficiently converged density at lambda=0.0_dp). 1187 else if( abs(lambda_new-lambda_predict)<abs(lambda_old-lambda_predict) & 1188 & .or. & 1189 & ( abs(lambda_old)<1.0d-6 .and. & 1190 & ilinmin>1 .and. & 1191 & etotal>etotal_previous ) & 1192 & )then 1193 f_fftgr(:,:,1)=vtrial(:,:) 1194 tmp=i_rhor(3) ; i_rhor(3)=i_rhor(2) ; i_rhor(2)=tmp 1195 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 1196 tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) 1197 i_vrespc(2)=i_vrespc(1); i_vrespc(1)=tmp; 1198 tmp=i_vresid(3); i_vresid(3)=i_vresid(2) 1199 i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 1200 if(moved_atm_inside==1)then 1201 f_atm(:,:,1)=xred(:,:) 1202 f_atm(:,:,i_rhor(2))=xred(:,:) 1203 end if 1204 d_lambda_old2=lambda_old-lambda_new 1205 lambda_old=lambda_new 1206 etotal_old=etotal 1207 resid_old=resid_new(1) 1208 d2edv2_old2=d2edv2_new 1209 dedv_old=dedv_new 1210 dedv_old2=dedv_new 1211 ! if(abs(lambda_new-lambda_predict)*2.0_dp& 1212 ! & /abs(lambda_new+lambda_predict) <reduction) end_linmin=1 1213 1214 ! Here lambda_old is the closest to lambda_predict (except for avoiding 1215 ! lambda_old==0.0_dp) 1216 else 1217 tmp=i_vresid(3) ; i_vresid(3)=i_vresid(1) ; i_vresid(1)=tmp 1218 f_fftgr(:,:,i_rhor(3))=rhor(:,:) 1219 if(moved_atm_inside==1) f_atm(:,:,i_rhor(3))=xred(:,:) 1220 tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(1) ; i_vrespc(1)=tmp 1221 d_lambda_old2=lambda_new-lambda_old 1222 etotal_previous=etotal 1223 d2edv2_old2=d2edv2_old 1224 dedv_old2=dedv_old 1225 ! if(abs(lambda_old-lambda_predict)*2.0_dp& 1226 ! & /abs(lambda_old+lambda_predict) <reduction) end_linmin=1 1227 end if 1228 1229 ! If the interval has not yet been sufficiently reduced, 1230 ! continue the search 1231 if(end_linmin==0)then 1232 lambda_new=lambda_predict 1233 1234 ! DEBUG 1235 ! write(std_out,'(a,2es16.6)' )& 1236 ! & ' scfcge : continue search, lambda_old,lambda_new=',lambda_old,lambda_new 1237 ! write(std_out,'(a,2es16.6)' )& 1238 ! & ' scfcge : f_fftgr(3:4,1,1)=',f_fftgr(3:4,1,1) 1239 ! write(std_out,'(a,2es16.6)' )& 1240 ! & ' scfcge : f_fftgr(3:4,1,6)=',f_fftgr(3:4,1,6) 1241 ! ENDDEBUG 1242 1243 vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6) 1244 if(moved_atm_inside==1)then 1245 xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6) 1246 end if 1247 1248 ilinmin=ilinmin+1 1249 ! 1250 ! Here generates a starting point for next line search 1251 else 1252 iline_cge=iline_cge+1 1253 if(end_linmin==1)ilinmin=0 1254 lambda_old=0.0_dp 1255 1256 ! In order to generate the new step, take into account previous 1257 ! optimal lambdas (including those of previous ion moves), 1258 ! and the selected new one, if it is positive. 1259 ! However, wait iline_cge>1 to select new ones. 1260 ! lambda_adapt has been initialized at 1.0_dp 1261 if(iline_cge>1 .and. lambda_new>0.0_dp )then 1262 ! Actually compute a geometric mean 1263 lambda_adapt= ( lambda_adapt**(dble(iline_cge-1)) * abs(lambda_new)) & 1264 & **(1.0_dp/dble(iline_cge)) 1265 ! In order to recover the previous algorithm, it is enough 1266 ! to decomment the next line 1267 ! lambda_adapt=1.0_dp 1268 end if 1269 lambda_new=lambda_adapt 1270 1271 vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2)) 1272 if(moved_atm_inside==1)then 1273 xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2)) 1274 end if 1275 1276 ! End choice between continue line minim and determine new direction 1277 end if 1278 1279 ! 1280 ! ------------------------------- 1281 1282 ! Here perform the CG step 1283 1284 else if(ilinmin==0)then 1285 1286 ! Compute the approximate energy derivatives dedv_mix,dedv_new,dedv_old 1287 choice=3 1288 call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,& 1289 & f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,& 1290 & natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred) 1291 1292 dedv_mix=dedv_mix/lambda_new 1293 dedv_new=dedv_new/lambda_new 1294 dedv_old=dedv_old/lambda_new 1295 1296 ! DEBUG 1297 ! write(message, '(a,3es12.4)' )' scfcge: lambda_adapt',& 1298 ! & lambda_adapt 1299 ! call wrtout(std_out,message,'COLL') 1300 1301 ! write(message, '(a,3es12.4)' )' scfcge: dedv_old,dedv_new,dedv_mix',& 1302 ! & dedv_old,dedv_new,dedv_mix 1303 ! call wrtout(std_out,message,'COLL') 1304 ! ENDDEBUG 1305 1306 ! Then, compute a predicted point, either along the line, 1307 ! or in a 2D plane 1308 testcg=1 1309 if(testcg==0)then 1310 ! This part corresponds to steepest descent, 1311 ! in which the line minimisation can be done 1312 ! using different algorithms, varying with the value of choice 1313 choice=1 1314 if(iscf==6)choice=2 1315 call findminscf(choice,dedv_new,dedv_old,dedv_predict,& 1316 & d2edv2_new,d2edv2_old,d2edv2_predict,& 1317 & etotal,etotal_old,etotal_predict,& 1318 & lambda_new,lambda_old,lambda_predict,errid_,message) 1319 if (errid_ /= AB7_NO_ERROR) then 1320 call wrtout(std_out,message,'COLL') 1321 end if 1322 lambda_predict2=0.0_dp 1323 ! Suppress the next line for debugging (there is another such line) 1324 status=0 1325 else 1326 ! This part corresponds to conjugate gradient 1327 ! A 2D minimisation is performed 1328 ! oldest direction is labelled 2 1329 ! newest direction is labelled 1 1330 de1=dedv_old ; de2=dedv_old2 1331 d2e11=(dedv_new-dedv_old)/lambda_new 1332 d2e22=d2edv2_old2 1333 d2e12=(dedv_mix-dedv_old)/d_lambda_old2 1334 ! The system to be solved is 1335 ! 0 = de1 + lambda1 d2e11 + lambda2 d2d12 1336 ! 0 = de2 + lambda1 d2e12 + lambda2 d2d22 1337 determ=d2e11*d2e22-d2e12*d2e12 1338 lambda_predict=-(de1*d2e22-de2*d2e12)/determ 1339 lambda_predict2=(de1*d2e12-de2*d2e11)/determ 1340 d2edv2_new=d2e11 ; d2edv2_old=d2e11 1341 end if 1342 1343 ! DEBUG 1344 ! write(message, '(a,5es11.3)' )' scfcge: de1,de2,d2e11,d2e22,d2e12',& 1345 ! & de1,de2,d2e11,d2e22,d2e12 1346 ! call wrtout(std_out,message,'COLL') 1347 ! write(std_out,'(a,2es12.4)' )' scfcge: la_predict,la_predict2',& 1348 ! & lambda_predict,lambda_predict2 1349 ! ----- 1350 ! write(std_out,*)'residues ', 1351 ! !$ de1+lambda_predict*d2e11+lambda_predict2*d2e12, 1352 ! !$ de2+lambda_predict*d2e12+lambda_predict2*d2e22 1353 ! if(.true.)stop 1354 ! ENDDEBUG 1355 ! 1356 1357 ! Determine the region of the 2D search space 1358 ! in which the predicted point is located, 1359 ! or use linear indicator to decide interpolation 1360 ! and advance to next 2D search. 1361 end_linmin=0 1362 write(message, '(a,2i3)' )' nlinear, ilinear',nlinear,ilinear 1363 call wrtout(std_out,message,'COLL') 1364 if(lambda_predict<0.0_dp)then 1365 ! Something is going wrong. Just take a reasonable step 1366 ! along the steepest descent direction (Region III). 1367 ! Actually, Region I and region III are treated in the same way later. 1368 ! In effect, this corresponds to restart the algorithm 1369 end_linmin=3 1370 ! Also put ilinear and nlinear to 0 1371 ilinear=0 1372 nlinear=0 1373 ! Decrease the adaptive step to predict next direction 1374 lambda_adapt=lambda_adapt*0.7_dp 1375 else if(ilinear>=1) then 1376 ! Region IV : will do an interpolation 1377 end_linmin=4 1378 ilinear=ilinear-1 1379 else if(abs(lambda_predict2)>reduction .or.& 1380 & lambda_predict<0.5_dp .or.& 1381 & lambda_predict>2.5_dp .or.& 1382 & lambda_predict-abs(lambda_predict2)/reduction <0.0_dp ) then 1383 ! Region II : lambda_predict is not too good, and not too bad. 1384 end_linmin=2 1385 else if (abs(1.0_dp-lambda_predict)<reduction)then 1386 ! Region I, the out-of-line point is OK. 1387 end_linmin=1 1388 else 1389 ! If everything fails, then region II. 1390 end_linmin=2 1391 end if 1392 1393 ! DEBUG 1394 ! write(message, '(a,2es12.4,i2)' )& 1395 ! & ' scfcge : la_predict, la_predict2, region',& 1396 ! & lambda_predict,lambda_predict2,end_linmin 1397 ! call wrtout(std_out,message,'COLL') 1398 ! ENDDEBUG 1399 1400 ! Treat region I, in the same way as region III 1401 if(end_linmin==1 .or. end_linmin==3)then 1402 1403 ! In region I, the line search is 1404 ! along vtrial-vtrial_old. 1405 ! The closest point is the new point 1406 ! thus to be transfered in the "old" locations 1407 1408 do isp=1,nspden 1409 do ifft=1,cplex*nfft 1410 f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new 1411 end do 1412 end do 1413 f_fftgr(:,:,1)=vtrial(:,:) 1414 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 1415 if(moved_atm_inside==1)then 1416 f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new 1417 f_atm(:,:,1)=xred(:,:) 1418 f_atm(:,:,i_rhor(2))=xred(:,:) 1419 end if 1420 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 1421 tmp=i_vresid(3) ; i_vresid(3)=i_vresid(2) 1422 i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 1423 d_lambda_old2=-lambda_new 1424 lambda_old=lambda_new 1425 etotal_old=etotal 1426 resid_old=resid_new(1) 1427 d2edv2_old=d2edv2_new 1428 dedv_old=dedv_new 1429 1430 ! Region I or III : one is close of the 2D minimum, 1431 ! or lambda_predict was negative (indicate a problem of convergence) 1432 ! Compute next trial potential along the 1433 ! PC residual and not along this search direction. 1434 ilinmin=0 1435 ! Question : isn t it here that one should prevent region I to called 1436 ! itself more than 1 time ??? 1437 ! Here the small difference between region I and region III 1438 if(end_linmin==3)ilinmin=1 1439 lambda_old=0.0_dp 1440 lambda_new=lambda_adapt 1441 1442 vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2)) 1443 if(moved_atm_inside==1)then 1444 xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2)) 1445 end if 1446 ! The new vtrial has been generated 1447 1448 else 1449 1450 ! Here region II or IV 1451 ilinmin=1 1452 if (lambda_predict==0._dp) then 1453 gamma=zero 1454 else 1455 gamma=lambda_predict2/lambda_predict 1456 end if 1457 ! Compute new search direction and trial potential 1458 write(message,*)' compute new search direction ' 1459 call wrtout(std_out,message,'COLL') 1460 do isp=1,nspden 1461 do ifft=1,cplex*nfft 1462 f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new+ & 1463 & gamma*f_fftgr(ifft,isp,6) 1464 end do 1465 end do 1466 vtrial(:,:)=f_fftgr(:,:,1)+ lambda_predict*f_fftgr(:,:,6) 1467 if(moved_atm_inside==1)then 1468 f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new+ gamma*f_atm(:,:,6) 1469 xred(:,:)=f_atm(:,:,1)+ lambda_predict*f_atm(:,:,6) 1470 end if 1471 1472 ! If end_linmin==2, then this vtrial is the good one 1473 1474 if(end_linmin==2)then 1475 1476 lambda_old=0.0_dp 1477 lambda_new=lambda_predict 1478 1479 else if(end_linmin==4)then 1480 1481 ! predict the result of the computation at the trial potential 1482 ! defined in the end_linmin==2 case 1483 gamma=lambda_predict2/d_lambda_old2 1484 ratio=lambda_predict/lambda_new 1485 1486 ! Take care of vtrial 1487 f_fftgr(:,:,1)=vtrial(:,:) 1488 1489 ABI_MALLOC(tmp_fft1,(cplex*nfft,nspden)) 1490 ! Take care of vresid 1491 tmp_fft1(:,:)=f_fftgr(:,:,i_vresid(2)) 1492 f_fftgr(:,:,i_vresid(2))=tmp_fft1(:,:)& 1493 & +ratio*(f_fftgr(:,:,i_vresid(1))-tmp_fft1(:,:))& 1494 & +gamma*(f_fftgr(:,:,i_vresid(3))-tmp_fft1(:,:)) 1495 f_fftgr(:,:,i_vresid(3))=tmp_fft1(:,:) 1496 1497 ! Take care of rhor 1498 tmp_fft1(:,:)=f_fftgr(:,:,i_rhor(2)) 1499 f_fftgr(:,:,i_rhor(2))=tmp_fft1(:,:)& 1500 & +ratio*(rhor(:,:)-tmp_fft1(:,:))& 1501 & +gamma*(f_fftgr(:,:,i_rhor(3))-tmp_fft1(:,:)) 1502 f_fftgr(:,:,i_rhor(3))=tmp_fft1(:,:) 1503 1504 ! Take care of vrespc 1505 tmp_fft1(:,:)=f_fftgr(:,:,i_vrespc(2)) 1506 f_fftgr(:,:,i_vrespc(2))=tmp_fft1(:,:)& 1507 & +ratio*(f_fftgr(:,:,i_vrespc(1))-tmp_fft1(:,:))& 1508 & +gamma*(f_fftgr(:,:,i_vrespc(3))-tmp_fft1(:,:)) 1509 f_fftgr(:,:,i_vrespc(3))=tmp_fft1(:,:) 1510 ABI_FREE(tmp_fft1) 1511 1512 if(moved_atm_inside==1)then 1513 do idir=1,3 1514 do iatom=1,natom 1515 1516 ! Take care of xred 1517 f_atm(idir,iatom,1)=xred(idir,iatom) 1518 1519 ! Take care of -HF forces 1520 temp=f_atm(idir,iatom,i_vresid(2)) 1521 f_atm(idir,iatom,i_vresid(2))=f_atm(idir,iatom,i_vresid(2))& 1522 & +ratio*(f_atm(idir,iatom,i_vresid(1))-f_atm(idir,iatom,i_vresid(2)))& 1523 & +gamma*(f_atm(idir,iatom,i_vresid(3))-f_atm(idir,iatom,i_vresid(2))) 1524 f_atm(idir,iatom,i_vresid(3))=temp 1525 1526 ! Take care of old xreds 1527 temp=f_atm(idir,iatom,i_rhor(2)) 1528 f_atm(idir,iatom,i_rhor(2))=f_atm(idir,iatom,i_rhor(2))& 1529 & +ratio*( xred(idir,iatom) -f_atm(idir,iatom,i_rhor(2)))& 1530 & +gamma*(f_atm(idir,iatom,i_rhor(3))-f_atm(idir,iatom,i_rhor(2))) 1531 f_atm(idir,iatom,i_rhor(3))=temp 1532 1533 ! Take care of preconditioned changes of atomic positions 1534 temp=f_atm(idir,iatom,i_vrespc(2)) 1535 f_atm(idir,iatom,i_vrespc(2))=f_atm(idir,iatom,i_vrespc(2))& 1536 & +ratio*(f_atm(idir,iatom,i_vrespc(1))-f_atm(idir,iatom,i_vrespc(2)))& 1537 & +gamma*(f_atm(idir,iatom,i_vrespc(3))-f_atm(idir,iatom,i_vrespc(2))) 1538 f_atm(idir,iatom,i_vrespc(3))=temp 1539 1540 end do 1541 end do 1542 end if 1543 1544 ! Since we are at the 2D minimum, the derivative is supposed 1545 ! to vanish. Note that dedv_old should not change, by contrast. 1546 dedv_old2=0.0_dp 1547 d_lambda_old2=-lambda_predict 1548 d2edv2_old2=-dedv_old/lambda_predict 1549 lambda_old=lambda_predict 1550 ilinmin=0 1551 1552 ! So, jump to the next line 1553 iline_cge=iline_cge+1 1554 write(message,*)' energy CG update : after 2D interpolation,' 1555 call wrtout(std_out,message,'COLL') 1556 write(message,*)' computation in the next plane ' 1557 call wrtout(std_out,message,'COLL') 1558 write(message,*) 1559 call wrtout(std_out,message,'COLL') 1560 lambda_old=0.0_dp 1561 lambda_new=lambda_adapt 1562 1563 vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2)) 1564 if(moved_atm_inside==1)then 1565 xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2)) 1566 end if 1567 1568 ! The new trial potential is now generated 1569 1570 ! End the specific treatment of region IV 1571 end if 1572 ! 1573 ! End the choice between treatment of region I, II, or IV 1574 end if 1575 1576 ! End of choice between initialisation or more developed parts of the CG algorithm 1577 else 1578 errid = AB7_ERROR_MIXING_ARG 1579 errmess = 'scfcge : BUG You should not be here ! ' 1580 return 1581 end if 1582 1583 !-------------------------------------- 1584 1585 !Write information : it will be easy to read by typing grep scfcge logfile 1586 1587 if(istep==1)then 1588 write(message,'(a,a,a)') ' scfcge:',ch10,' scfcge:istep-iline_cge-ilinmin lambda etot resid ' 1589 call wrtout(std_out,message,'COLL') 1590 end if 1591 1592 if(ilinmin_input/=0 .or. istep==1)then 1593 ! Usual line minimisation step 1594 1595 if(iline_cge_input<10)then 1596 write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' )& 1597 & ' scfcge: actual ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input 1598 else 1599 write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' )& 1600 & ' scfcge: actual ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input 1601 end if 1602 call wrtout(std_out,message,'COLL') 1603 1604 if( (end_linmin==1.or.end_linmin==-1) .and. istep/=1 )then 1605 1606 if(end_linmin==1)then 1607 write(message, '(a,es13.4,a,i2,a,a)' )& 1608 & ' scfcge: predict ',lambda_predict,& 1609 & ' suff. close => next line, ilinear=',ilinear,ch10,& 1610 & ' scfcge:' 1611 else if(end_linmin==-1)then 1612 write(message, '(a,es13.4,a,a,a)' )& 1613 & ' scfcge: predict ',lambda_predict,& 1614 & ' restart the algorithm ',ch10,& 1615 & ' scfcge:' 1616 end if 1617 call wrtout(std_out,message,'COLL') 1618 1619 if(iline_cge_input<9)then 1620 write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' ) & 1621 & ' scfcge: start ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old 1622 else 1623 write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' ) & 1624 & ' scfcge: start ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old 1625 end if 1626 call wrtout(std_out,message,'COLL') 1627 1628 else if(istep/=1) then 1629 write(message, '(a,es13.4,a)' )& 1630 & ' scfcge: predict ',lambda_predict,& 1631 & ' not close enough => continue minim.' 1632 call wrtout(std_out,message,'COLL') 1633 end if 1634 1635 else 1636 ! CG prediction 1637 if(iline_cge_input<10)then 1638 write(message, '(a,i4,a,i1,a,es11.4,es20.12,es12.4,a,i1)' )& 1639 & ' scfcge: actual ',istep,'-',iline_cge_input,'-off',& 1640 & lambda_adapt,etotal_input,resid_input,', end=',end_linmin 1641 else 1642 write(message, '(a,i3,a,i2,a,es11.4,es20.12,es12.4,a,i1)' )& 1643 & ' scfcge: actual ',istep,'-',iline_cge_input,'-off',& 1644 & lambda_adapt,etotal_input,resid_input,', end=',end_linmin 1645 end if 1646 call wrtout(std_out,message,'COLL') 1647 1648 if(end_linmin==4)then 1649 write(message, '(a)' ) ' scfcge:' 1650 call wrtout(std_out,message,'COLL') 1651 end if 1652 1653 end if 1654 1655 end subroutine scfcge
m_ab7_mixing/scfopt [ Functions ]
[ Top ] [ m_ab7_mixing ] [ Functions ]
NAME
scfopt
FUNCTION
Compute the next vtrial of the SCF cycle. Possible algorithms are : simple mixing, Anderson (order 1 or 2), Pulay
INPUTS
cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX iscf= 2 => simple mixing = 3,4 => Anderson mixing = 7 => Pulay mixing istep= number of the step in the SCF cycle mpicomm=the mpi communicator used for the summation comm_atom=the mpi communicator over atoms ; PAW only (optional argument) mpi_summarize=set it to .true. if parallelisation is done over FFT nfft=(effective) number of FFT grid points (for this processor) npawmix=-PAW only- number of spherical part elements to be mixed nspden=number of spin-density components n_fftgr=third dimension of the array f_fftgr n_index=dimension for indices of potential/density (see ivrespc, i_vtrial...) opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential 1 vtrial (and also f_fftgr) actually contains the trial density pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part usepaw= 0 for non paw calculation; =1 for paw calculation
OUTPUT
(see side effects)
SIDE EFFECTS
vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
the input preconditioned residual potential
at output, it is the new trial potential .
f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
The input vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(1)).
The old vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(2)).
The input preconditioned residual potential is in f_fftgr(:,:,i_vrespc(1))
Two input old preconditioned residual potentials in f_fftgr(:,:,i_vrespc(2)) and f_fftgr(:,:,i_vrespc(3))
Before output a permutation of i_vrespc(1), i_vrespc(2) and i_vrespc(3) occurs, without
actually copying all the data (change of pointer).
i_vrespc(n_index)=index of the preconditioned residual potentials (present and past) in the array f_fftgr
i_vtrial(n_index) =indices of the potential (present and past) in the array f_fftgr
==== if usepaw==1
f_paw(npawmix,n_fftgr*mffmem*usepaw)=different functions used for PAW
(same as f_fftgr but for spherical part)
vpaw(npawmix*usepaw)=at input, the aug. occupancies (rhoij) that gave
the input preconditioned residual potential
at output, it is the new aug. occupancies.
SOURCE
1868 subroutine scfopt(cplex,f_fftgr,f_paw,iscf,istep,i_vrespc,i_vtrial,& 1869 & mpicomm,mpi_summarize,nfft,npawmix,nspden,n_fftgr,& 1870 & n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial,errid,errmess, & 1871 & comm_atom) ! optional 1872 1873 !Arguments ------------------------------------ 1874 !scalars 1875 integer,intent(in) :: cplex,iscf,istep,n_fftgr,n_index,nfft 1876 integer,intent(in) :: npawmix,nspden,opt_denpot,pawoptmix,usepaw,mpicomm 1877 integer, intent(in),optional :: comm_atom 1878 integer,intent(out) :: errid 1879 character(len = 500), intent(out) :: errmess 1880 logical, intent(in) :: mpi_summarize 1881 real(dp), intent(out) :: vresid 1882 !arrays 1883 integer,intent(inout) :: i_vrespc(n_index),i_vtrial(n_index) 1884 real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr) 1885 real(dp),intent(inout) :: f_paw(npawmix,n_fftgr*usepaw),vpaw(npawmix*usepaw) 1886 real(dp),intent(inout) :: vtrial(cplex*nfft,nspden) 1887 !Local variables------------------------------- 1888 !scalars 1889 integer,parameter :: npulaymax=50 1890 integer :: i_vstore,ierr,ifft,ii,index,isp,jj,comm_atom_,niter,npulay,tmp 1891 real(dp),save :: prod_resid_old,resid_old,resid_old2 1892 real(dp) :: aa1,aa2,bb,cc1,cc2,current,det,lambda,lambda2,resid_best 1893 character(len=500) :: message 1894 !arrays 1895 integer,allocatable :: ipiv(:) 1896 real(dp),save :: amat(npulaymax+1,npulaymax+1) 1897 real(dp) :: mpibuff(2),prod_resid(1),prod_resid2(1),resid_new(1) 1898 real(dp),allocatable :: alpha(:),amatinv(:,:),amat_paw(:),rwork(:) 1899 1900 ! ************************************************************************* 1901 1902 !DEBUG 1903 !write(std_out,*)' scfopt : enter ; istep,iscf ',istep,iscf 1904 !ENDDEBUG 1905 1906 errid = AB7_NO_ERROR 1907 1908 comm_atom_=xmpi_comm_self; if(present(comm_atom)) comm_atom_=comm_atom 1909 1910 i_vstore=i_vtrial(1) 1911 if (iscf==4) i_vstore=i_vtrial(2) 1912 if (iscf==7) then 1913 if (modulo(n_fftgr, 2) == 0 ) then 1914 npulay=(n_fftgr-2)/2 1915 else 1916 npulay=(n_fftgr-1)/2 1917 end if 1918 i_vstore=i_vtrial(npulay) 1919 else 1920 npulay=0 1921 end if 1922 1923 !Compute the new residual resid_new, from f_fftgr/f_paw(:,:,i_vrespc(1)) 1924 call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr) 1925 if (usepaw==1.and.pawoptmix==1) then 1926 do index=1,npawmix 1927 resid_new(1)=resid_new(1)+f_paw(index,i_vrespc(1))**2 1928 end do 1929 call xmpi_sum(resid_new(1),comm_atom_,ierr) 1930 end if 1931 vresid = resid_new(1) 1932 1933 !_______________________________________________________________ 1934 !Here use only the preconditioning, or initialize the other algorithms 1935 1936 if (istep==1 .or. iscf==2) then 1937 write(message,'(2a)') ch10,' Simple mixing update:' 1938 call wrtout(std_out,message,'COLL') 1939 1940 write(message,*)' residual square of the potential: ',resid_new(1) 1941 call wrtout(std_out,message,'COLL') 1942 1943 ! Store information for later use 1944 if (iscf==3.or.iscf==4) resid_old=resid_new(1) 1945 if (iscf==7) then 1946 amat(:,:)=zero 1947 amat(1,1)=resid_new(1) 1948 end if 1949 1950 ! Compute new vtrial (and new rhoij if PAW) 1951 if (iscf/=2) f_fftgr(:,:,i_vstore)=vtrial(:,:) 1952 vtrial(:,:)=vtrial(:,:)+f_fftgr(:,:,i_vrespc(1)) 1953 if (usepaw==1) then 1954 if (iscf/=2) f_paw(:,i_vstore)=vpaw(:) 1955 vpaw(:)=vpaw(:)+f_paw(:,i_vrespc(1)) 1956 end if 1957 1958 ! _______________________________________________________________ 1959 ! Here Anderson algorithm using one previous iteration 1960 else if((istep==2 .or. iscf==3).and.iscf/=7)then 1961 1962 write(message,'(2a)') ch10,' Anderson update:' 1963 call wrtout(std_out,message,'COLL') 1964 1965 write(message,*)' residual square of the potential: ',resid_new(1) 1966 call wrtout(std_out,message,'COLL') 1967 1968 ! Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2)) 1969 call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,& 1970 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 1971 if (usepaw==1.and.pawoptmix==1) then 1972 do index=1,npawmix 1973 prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2)) 1974 end do 1975 call xmpi_sum(prod_resid(1),comm_atom_,ierr) 1976 end if 1977 1978 ! Compute mixing factor 1979 lambda=(resid_new(1)-prod_resid(1))/(resid_new(1)+resid_old-2*prod_resid(1)) 1980 write(message,*)' mixing of old trial potential: ',lambda 1981 call wrtout(std_out,message,'COLL') 1982 1983 ! Evaluate best residual square on the line 1984 resid_best=(1.0_dp-lambda)*(1.0_dp-lambda)*resid_new(1)& 1985 & +(1.0_dp-lambda)*lambda *2*prod_resid(1)& 1986 & +lambda *lambda *resid_old 1987 write(message,*)' predicted best residual square on the line: ',resid_best 1988 call wrtout(std_out,message,'COLL') 1989 1990 ! Store information for later use 1991 if (iscf==4) then 1992 prod_resid_old=prod_resid(1) 1993 resid_old2=resid_old 1994 end if 1995 resid_old=resid_new(1) 1996 1997 ! Save latest trial potential and compute new trial potential 1998 do isp=1,nspden 1999 do ifft=1,cplex*nfft 2000 current=vtrial(ifft,isp) 2001 vtrial(ifft,isp)=(one-lambda)*(current +f_fftgr(ifft,isp,i_vrespc(1)))& 2002 & +lambda *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2))) 2003 f_fftgr(ifft,isp,i_vstore)=current 2004 end do 2005 end do 2006 2007 ! PAW: save latest rhoij and compute new rhoij 2008 do index=1,npawmix 2009 current=vpaw(index) 2010 vpaw(index)=(one-lambda)*(current +f_paw(index,i_vrespc(1)))& 2011 & +lambda *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2))) 2012 f_paw(index,i_vstore)=current 2013 end do 2014 2015 ! _______________________________________________________________ 2016 ! Here Anderson algorithm using two previous iterations 2017 else if(iscf==4.and.iscf/=7)then 2018 2019 write(message,'(2a)') ch10,' Anderson (order 2) update:' 2020 call wrtout(std_out,message,'COLL') 2021 2022 write(message,*)' residual square of the potential: ',resid_new(1) 2023 call wrtout(std_out,message,'COLL') 2024 2025 ! Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2)) 2026 call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,& 2027 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 2028 if (usepaw==1.and.pawoptmix==1) then 2029 do index=1,npawmix 2030 prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2)) 2031 end do 2032 end if 2033 2034 ! Compute prod_resid2 from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(3)) 2035 call dotprodm_v(cplex,1,prod_resid2,i_vrespc(1),i_vrespc(3),mpicomm,mpi_summarize,1,1,& 2036 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 2037 if (usepaw==1.and.pawoptmix==1) then 2038 do index=1,npawmix 2039 prod_resid2(1)=prod_resid2(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(3)) 2040 end do 2041 ! MPI reduction 2042 mpibuff(1)=prod_resid(1);mpibuff(2)=prod_resid2(1) 2043 call xmpi_sum(mpibuff,comm_atom_,ierr) 2044 prod_resid(1)=mpibuff(1);prod_resid2(1)=mpibuff(2) 2045 end if 2046 2047 ! Compute mixing factors 2048 aa1=resid_new(1)+resid_old -two*prod_resid (1) 2049 aa2=resid_new(1)+resid_old2-two*prod_resid2(1) 2050 bb =resid_new(1)+prod_resid_old-prod_resid(1)-prod_resid2(1) 2051 cc1=resid_new(1)-prod_resid (1) 2052 cc2=resid_new(1)-prod_resid2(1) 2053 det=aa1*aa2-bb*bb 2054 lambda =(aa2*cc1-bb*cc2)/det 2055 lambda2=(aa1*cc2-bb*cc1)/det 2056 write(message,*)' mixing of old trial potentials: ',lambda,lambda2 2057 call wrtout(std_out,message,'COLL') 2058 2059 ! Store information for later use 2060 prod_resid_old=prod_resid(1) 2061 resid_old2=resid_old 2062 resid_old=resid_new(1) 2063 2064 ! Save latest trial potential and compute new trial potential 2065 do isp=1,nspden 2066 do ifft=1,cplex*nfft 2067 current=vtrial(ifft,isp) 2068 vtrial(ifft,isp)=& 2069 & (one-lambda-lambda2)*(current +f_fftgr(ifft,isp,i_vrespc(1)))& 2070 & +lambda *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))& 2071 & +lambda2 *(f_fftgr(ifft,isp,i_vtrial(2))+f_fftgr(ifft,isp,i_vrespc(3))) 2072 f_fftgr(ifft,isp,i_vstore)=current 2073 end do 2074 end do 2075 2076 ! PAW: save latest rhoij and compute new rhoij 2077 do index=1,npawmix 2078 current=vpaw(index) 2079 vpaw(index)=& 2080 & (one-lambda-lambda2)*(current +f_paw(index,i_vrespc(1)))& 2081 & +lambda *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))& 2082 & +lambda2 *(f_paw(index,i_vtrial(2))+f_paw(index,i_vrespc(3))) 2083 f_paw(index,i_vstore)=current 2084 end do 2085 2086 ! _______________________________________________________________ 2087 ! Here Pulay algorithm 2088 else if(iscf==7)then 2089 2090 niter=min(istep,npulay+1) 2091 2092 write(message,'(2a,i2,a)') ch10,' Pulay update with ',niter-1,' previous iterations:' 2093 call wrtout(std_out,message,'COLL') 2094 2095 if (npulay>npulaymax) then 2096 errid = AB7_ERROR_MIXING_CONVERGENCE 2097 write(errmess, '(4a)' ) ch10,& 2098 & ' scfopt: ERROR - ',ch10,& 2099 & ' Too many iterations required for Pulay algorithm (<50) !' 2100 return 2101 end if 2102 2103 ! Compute "A" matrix 2104 if (istep>npulay+1) then 2105 do jj=1,niter-1 2106 do ii=1,niter-1 2107 amat(ii,jj)=amat(ii+1,jj+1) 2108 end do 2109 end do 2110 end if 2111 if (usepaw==1.and.pawoptmix==1) then 2112 ABI_MALLOC(amat_paw,(niter)) 2113 amat_paw(:)=zero 2114 do ii=1,niter 2115 do index=1,npawmix 2116 amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(1+niter-ii)) 2117 end do 2118 end do 2119 call xmpi_sum(amat_paw,comm_atom_,ierr) 2120 end if 2121 do ii=1,niter 2122 call dotprodm_v(cplex,1,amat(ii,niter),i_vrespc(1),i_vrespc(1+niter-ii),mpicomm,mpi_summarize,1,1,& 2123 & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr) 2124 if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii) 2125 if (ii<niter) amat(niter,ii)=amat(ii,niter) 2126 end do 2127 if (usepaw==1.and.pawoptmix==1)then 2128 ABI_FREE(amat_paw) 2129 end if 2130 2131 ! Invert "A" matrix 2132 ABI_MALLOC(amatinv,(niter,niter)) 2133 amatinv(1:niter,1:niter)=amat(1:niter,1:niter) 2134 ABI_MALLOC(ipiv,(niter)) 2135 ABI_MALLOC(rwork,(niter)) 2136 call dgetrf(niter,niter,amatinv,niter,ipiv,ierr) 2137 call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr) 2138 ABI_FREE(ipiv) 2139 ABI_FREE(rwork) 2140 2141 ! Compute "alpha" factors 2142 ABI_MALLOC(alpha,(niter)) 2143 alpha=zero 2144 det=zero 2145 do ii=1,niter 2146 do jj=1,niter 2147 alpha(ii)=alpha(ii)+amatinv(jj,ii) 2148 det=det+amatinv(jj,ii) 2149 end do 2150 end do 2151 alpha(:)=alpha(:)/det 2152 ABI_FREE(amatinv) 2153 write(message,'(a,5(1x,g10.3))')' mixing of old trial potential: alpha(m:m-4)=',(alpha(ii),ii=niter,max(1,niter-4),-1) 2154 call wrtout(std_out,message,'COLL') 2155 2156 ! Save latest trial potential and compute new trial potential 2157 do isp=1,nspden 2158 do ifft=1,cplex*nfft 2159 current=vtrial(ifft,isp) 2160 vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,i_vrespc(1))) 2161 do ii=niter-1,1,-1 2162 vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) & 2163 & *(f_fftgr(ifft,isp,i_vtrial(niter-ii))+f_fftgr(ifft,isp,i_vrespc(1+niter-ii))) 2164 end do 2165 f_fftgr(ifft,isp,i_vstore)=current 2166 end do 2167 end do 2168 2169 ! PAW: save latest rhoij and compute new rhoij 2170 do index=1,npawmix 2171 current=vpaw(index) 2172 vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1))) 2173 do ii=niter-1,1,-1 2174 vpaw(index)=vpaw(index)+alpha(ii) & 2175 & *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii))) 2176 end do 2177 f_paw(index,i_vstore)=current 2178 end do 2179 2180 ABI_FREE(alpha) 2181 ! _______________________________________________________________ 2182 ! End of choice of optimization method 2183 end if 2184 2185 !Permute potential indices 2186 if (iscf==3) then 2187 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 2188 else if (iscf==4) then 2189 tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 2190 tmp=i_vtrial(2) ; i_vtrial(2)=i_vtrial(1) ; i_vtrial(1)=tmp 2191 else if (iscf==7) then 2192 tmp=i_vtrial( npulay) 2193 do ii= npulay,2,-1 2194 i_vtrial(ii)=i_vtrial(ii-1) 2195 end do 2196 i_vtrial(1)=tmp 2197 tmp=i_vrespc(1+npulay) 2198 do ii=1+npulay,2,-1 2199 i_vrespc(ii)=i_vrespc(ii-1) 2200 end do 2201 i_vrespc(1)=tmp 2202 end if 2203 2204 end subroutine scfopt