From 5bb7336871e9b551e5b7d1f20fac5c4725873d62 Mon Sep 17 00:00:00 2001 From: Panadestein Date: Tue, 19 Jan 2021 16:33:56 +0100 Subject: [PATCH] Added CHAMP jastrow caller --- hpsi.f | 229 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 229 insertions(+) create mode 100644 hpsi.f diff --git a/hpsi.f b/hpsi.f new file mode 100644 index 0000000..e8ba28c --- /dev/null +++ b/hpsi.f @@ -0,0 +1,229 @@ + subroutine hpsi(coord,psid,psij,energy,ipass,ifr) + +c Written by Cyrus Umrigar, modified by Claudia Filippi and A. Scemama +c modified by Claudio Amovilli and Franca Floris for PCM and QM-MMPOl + + implicit real*8(a-h,o-z) + include 'vmc.h' + include 'pseudo.h' + include 'force.h' + include 'pcm.h' + include 'mmpol.h' + include 'efield.h' + include 'optjas.h' + include 'mstates.h' + +c Calculates energy + + common /const/ pi,hb,etrial,delta,deltai,fbias,nelec,imetro,ipr + common /contr2/ ijas,icusp,icusp2,isc,ianalyt_lap + &,ifock,i3body,irewgt,iaver,istrch + common /elec/ nup,ndn + + common /coefs/ coef(MBASIS,MORB,MWF),nbasis,norb + + common /dets/ cdet(MDET,MSTATES,MWF),ndet + common /csfs/ ccsf(MDET,MSTATES,MWF),cxdet(MDET*MDETCSFX) + &,icxdet(MDET*MDETCSFX),iadet(MDET),ibdet(MDET),ncsf,nstates + + common /optwf_contrl/ ioptjas,ioptorb,ioptci,nparm + + common /wfsec/ iwftype(MFORCE),iwf,nwftype + + common /distance/ rshift(3,MELEC,MCENT),rvec_en(3,MELEC,MCENT) + &,r_en(MELEC,MCENT),rvec_ee(3,MMAT_DIM2),r_ee(MMAT_DIM2) + + common /slater/ slmui(MMAT_DIM),slmdi(MMAT_DIM) + &,fpu(3,MMAT_DIM),fpd(3,MMAT_DIM) + &,fppu(MMAT_DIM),fppd(MMAT_DIM) + &,ddx(3,MELEC),d2dx2(MELEC) + common /multislater/ detu(MDET),detd(MDET) + + common /jaso/ fso(MELEC,MELEC),fijo(3,MELEC,MELEC) + &,d2ijo(MELEC,MELEC),d2o,fsumo,fjo(3,MELEC) + + common /pseudo/ vps(MELEC,MCENT,MPS_L),vpso(MELEC,MCENT,MPS_L,MFORCE) + &,lpot(MCTYPE),nloc + + common /casula/ t_vpsp(MCENT,MPS_QUAD,MELEC),icasula,i_vpsp + common /pcm_hpsi/ pepcms,pepcmv,qopcm,enfpcm(MCHS) + common /mmpol_hpsi/ peQMdp,peQMq,eek_pol(3,MCHMM) + + common /velocity_jastrow/vj(3,MELEC),vjn(3,MELEC) + + common /multidet/ kref,numrep_det(MDET,2),irepcol_det(MELEC,MDET,2),ireporb_det(MELEC,MDET,2) + & ,iwundet(MDET,2),iactv(2),ivirt(2) + + common /Bloc/ b(MORB,MELEC),xmatu(MELEC**2),xmatd(MELEC**2) + & ,tildem(MELEC,MORB,2) + + common /ycompact/ ymat(MORB,MELEC,2,MSTATES),dymat(MORB,MELEC,2,MSTATES) + + common /force_analy/ iforce_analy + + ! ====RLPB for debug=== + common /config/ xold(3, MELEC), xnew(3, MELEC), vold(3, MELEC) + dimension ran(3, MELEC) + ! ===================== + + dimension coord(3,*),psid(*),energy(*) + dimension denergy(MSTATES),eloc_det(MDET,2),vpsp_det(2),dvpsp_dj(MPARMJ) + + iwf=iwftype(ifr) + +c pe_ee computed in distances (called from distances if iperiodic != 0) +c pe_en(loc) computed in distances if nloc=0 or called from distances if iperiodic!=0 +c pe_en(loc) computed in nonloc_pot if nloc!=0 and iperiodic=0 +c pe_en(nonloc) computed in nonloc_pot if nloc !=0 + +c distances needed for Jastrow, determinants, and potential energy + call distances(0,coord) + +c local potential contributions + call pot_local(pe_local) + +c external potential on a grid (e.g. MM from CPMD) + call p2gtid('qmmm:iqmmm',iqmmm,0,1) + if(iqmmm.eq.1) then + ext_pot=0 + call qmmm_extpot_ene(coord,nelec,ext_pot) + pe_local=pe_local+ext_pot + endif + +c external charges + if(iefield.eq.1) then + ext_pot=0 + call efield_extpot_ene(coord,nelec,ext_pot) + pe_local=pe_local+ext_pot + endif + +c PCM polarization charges + if(ipcm.gt.1) then + pepcms=0 + pepcmv=0 + call pcm_extpot_ene(coord,nelec,pepcms,pepcmv) + pepcm=pepcms+pepcmv + pe_local=pe_local+pepcm + endif + +c QM-MMPOL (charges+induced dipoles) + if(immpol.gt.1) then + peQMdp=0 + peQMq=0 + call mmpol_extpot_ene(coord,nelec,peQMdp,peQMq) + peQM=peQMdp+peQMq + pe_local=pe_local+peQM + endif + + if(ipr.ge.3) write(6,'(''pe_loc before nonloc_pot'',9f12.5)') pe_local + +c get contribution from jastrow (also compute derivatives wrt parameters and nuclei) + if(ianalyt_lap.eq.1) then + ! BEGIN DEBUG JASTROW RLPB + ! Here we are messing up with the xold variable in common block! + + ! We compare the output of both procedures + !call random_seed() + !do i = 1, 10 + ! call random_number(ran) + ! xold = ran + ! call distances(0, xold) + ! call jastrow(xold,vj,d2j,psij,ifr) + ! call irp_full_jastrow + !end do + + ! Here we benchmark both procedures + + call cpu_time(tchamp_init) + do i = 1, 100000 + call distances(0, coord) + call jastrow(coord,vj,d2j,psij,ifr) + end do + call cpu_time(tchamp_final) + + call cpu_time(tirp_init) + do i = 1, 100000 + call irp_full_jastrow + end do + call cpu_time(tirp_final) + + print *, "BEGIN ELAPSED TIMES" + print *, "CHAMP = ", tchamp_final - tchamp_init + print *, "IRP = ", tirp_final - tirp_init + print *, "END ELAPSED TIMES" + ! END DEBUG JASTROW RLPB + + !call jastrow(coord,vj,d2j,psij,ifr) + else + call jastrow_num(coord,vj,d2j,psij) + endif + if(ipr.ge.3) write(6,'(''d2j,psij'',9f12.5)') d2j,psij + +c compute reference determinant, its derivatives, and kinetic contribution to B_eloc and its derivatives + call determinant(ipass,coord,rvec_en,r_en) + + call compute_bmatrices_kin + +c compute pseudo-potential contribution +c nonloc_pot must be called after determinant because slater matrices are needed + if(nloc.gt.0) + & call nonloc_pot(coord,rshift,rvec_en,r_en,pe_local,vpsp_det,dvpsp_dj,t_vpsp,i_vpsp,ifr) + + if(ipr.ge.3) then + write(6,'(''pe_loc after nonloc_pot'',9f12.5)') pe_local + write(6,'(''pe_ref after nonloc_pot'',9f12.5)') (vpsp_det(ii),ii=1,2) + endif + + call multideterminant_hpsi(vj,vpsp_det,eloc_det) + + e_other=pe_local-hb*d2j + do 10 i=1,nelec + 10 e_other=e_other-hb*(vj(1,i)**2+vj(2,i)**2+vj(3,i)**2) + + do 30 istate=1,nstates + +c combine determinantal quantities to obtain trial wave function + call determinant_psit(psid(istate),istate) + +c compute energy using Ymat + denergy(istate)=0 + do 20 iab=1,2 + nel=nup + if(iab.eq.2) nel=ndn + do 20 jrep=ivirt(iab),norb + do 20 irep=iactv(iab),nel + 20 denergy(istate)=denergy(istate)+ymat(jrep,irep,iab,istate)*tildem(irep,jrep,iab) + denergy(istate)=denergy(istate)*detu(kref)*detd(kref)/psid(istate) + + energy(istate)=denergy(istate)+eloc_det(kref,1)+eloc_det(kref,2)+e_other + + if(ipr.ge.2) then + write(6,'(''state'',i4)') istate + write(6,'(''psid,psij'',9d12.5)') psid(istate),psij + write(6,'(''psitot '',e18.11)') psid(istate)*exp(psij) +c do k=1,ndet +c write(6,'(''psitot_k '',i6,3e18.8)') k, detu(k),detd(k),detu(k)*detd(k)*exp(psij) +c write(6,'(''psitot_k '',i6,3e18.8)') k, detu(k),detd(k),cdet(k,1,1)*detu(k)*detd(k)*exp(psij) +c enddo +c do 25 i=1,nelec +c do 25 k=1,3 +c 25 write(6,'(''vj'',2e18.11)') vj(k,i) + if(ipr.ge.3) write(6,'(''energy'',9f16.10)') energy(istate) + endif + + 30 continue + + if(ifr.eq.1) then + if(iforce_analy.eq.1) call compute_force(psid(1),denergy(1)) + + call optorb_compute(psid,energy,denergy) + + call optjas_deloc(psid,energy,dvpsp_dj,vj) + + call optci_deloc(eloc_det,e_other,psid,energy) + + call prop_compute(coord) + endif + + return + end