mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-29 06:35:39 -04:00
New implementation of QM/MM.
....EJB
This commit is contained in:
Eric Bylaska
parent
5a1523c213
commit
348ab26487
35 changed files with 1500 additions and 2899 deletions
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: c_cgsd_energy_gradient.F,v 1.4 2004-11-29 16:05:17 bylaska Exp $ c
|
||||
c $Id: c_cgsd_energy_gradient.F,v 1.5 2007-03-22 20:46:19 bylaska Exp $ c
|
||||
|
||||
* *******************************
|
||||
* * *
|
||||
|
|
@ -29,7 +29,7 @@ c $Id: c_cgsd_energy_gradient.F,v 1.4 2004-11-29 16:05:17 bylaska Exp $
|
|||
real*8 ion_rion
|
||||
external cpsp_semicore
|
||||
external control_allow_translation,ion_q_FixIon
|
||||
external ion_aname,pspw_qmmm_aname
|
||||
external ion_aname
|
||||
external ion_katm,ion_nion,control_version
|
||||
external ion_rion
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: ion.F,v 1.45 2007-03-19 19:05:35 bylaska Exp $
|
||||
* $Id: ion.F,v 1.46 2007-03-22 20:46:19 bylaska Exp $
|
||||
*
|
||||
|
||||
* ***************************
|
||||
|
|
@ -27,6 +27,7 @@
|
|||
value = value.and.MA_free_heap(atom_qm(2))
|
||||
|
||||
value = value.and.MA_free_heap(dti(2))
|
||||
value = value.and.MA_free_heap(amass(2))
|
||||
if (.not. value) call errquit('ion_end:free heap ',0, MA_ERR)
|
||||
|
||||
return
|
||||
|
|
@ -196,9 +197,9 @@
|
|||
|
||||
logical control_frac_coord,control_init_velocities
|
||||
logical parseqmmm
|
||||
double precision control_ion_time_step,ion_amass
|
||||
double precision control_ion_time_step,ion_amass_geom
|
||||
external control_frac_coord,control_init_velocities
|
||||
external control_ion_time_step,ion_amass
|
||||
external control_ion_time_step,ion_amass_geom
|
||||
external parseqmmm
|
||||
|
||||
call Parallel_taskid(taskid)
|
||||
|
|
@ -231,6 +232,8 @@
|
|||
> MA_alloc_get(mt_dbl,(3*nion),'r0',r0(2),r0(1))
|
||||
value = value.and.
|
||||
> MA_alloc_get(mt_dbl,(nion),'dti',dti(2),dti(1))
|
||||
value = value.and.
|
||||
> MA_alloc_get(mt_dbl,(nion),'amass',amass(2),amass(1))
|
||||
if (.not. value) call errquit('out of heap memory',0, MA_ERR)
|
||||
|
||||
call dcopy(3*nion,0.0d0,0,dbl_mb(r2(1)),1)
|
||||
|
|
@ -269,7 +272,8 @@
|
|||
* **** define dti ****
|
||||
dt = control_ion_time_step()
|
||||
do i=1,nion
|
||||
dbl_mb(dti(1)+i-1) = dt*dt/(ion_amass(i))
|
||||
dbl_mb(dti(1) +i-1) = dt*dt/(ion_amass_geom(i))
|
||||
dbl_mb(amass(1)+i-1) = ion_amass_geom(i)
|
||||
end do
|
||||
|
||||
c value = geom_destroy(geom)
|
||||
|
|
@ -415,10 +419,10 @@ c eki1=0.5d0*eki1
|
|||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_amass *
|
||||
* * ion_amass_geom *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function ion_amass(i)
|
||||
real*8 function ion_amass_geom(i)
|
||||
implicit none
|
||||
integer i
|
||||
|
||||
|
|
@ -436,10 +440,46 @@ c eki1=0.5d0*eki1
|
|||
if(.not.geom_mass_get(geom,i,mass))
|
||||
> call errquit(' geom_mass_get failed ',i, GEOM_ERR)
|
||||
|
||||
ion_amass = mass*1822.89d0
|
||||
ion_amass_geom = mass*1822.89d0
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_amass *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function ion_amass(i)
|
||||
implicit none
|
||||
integer i
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "errquit.fh"
|
||||
#include "ion.fh"
|
||||
|
||||
ion_amass = dbl_mb(amass(1)+i-1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_amass_ptr *
|
||||
* * *
|
||||
* ***************************
|
||||
integer function ion_amass_ptr()
|
||||
implicit none
|
||||
integer i
|
||||
|
||||
#include "ion.fh"
|
||||
|
||||
ion_amass_ptr = amass(1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_rcovalent *
|
||||
|
|
@ -528,11 +568,11 @@ c eki1=0.5d0*eki1
|
|||
* ***********************
|
||||
real*8 function ion_q(i)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
integer i
|
||||
|
||||
#include "stdio.fh"
|
||||
#include "geom.fh"
|
||||
#include "errquit.fh"
|
||||
|
||||
|
||||
integer geom
|
||||
|
|
@ -1118,6 +1158,23 @@ c eki1=0.5d0*eki1
|
|||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_katm_ptr *
|
||||
* * *
|
||||
* ***************************
|
||||
|
||||
integer function ion_katm_ptr()
|
||||
implicit none
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "ion.fh"
|
||||
|
||||
ion_katm_ptr = katm(1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_katm_qm *
|
||||
|
|
@ -1282,6 +1339,24 @@ c eki1=0.5d0*eki1
|
|||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_rion_ptr *
|
||||
* * *
|
||||
* ***************************
|
||||
|
||||
integer function ion_rion_ptr()
|
||||
implicit none
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "ion.fh"
|
||||
|
||||
ion_rion_ptr = r1(1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * ion_rion2 *
|
||||
|
|
|
|||
|
|
@ -3,11 +3,12 @@
|
|||
logical sym_inversion
|
||||
integer nion_mm,nion_qm
|
||||
integer nion,nkatm,nkatm_qm
|
||||
integer r2(2),r1(2),r0(2)
|
||||
integer r2(2),r1(2),r0(2),amass(2)
|
||||
integer katm(2),natm(2),atom(2)
|
||||
integer katm_qm(2),natm_qm(2),atom_qm(2)
|
||||
integer dti(2)
|
||||
common / ION / r2,r1,r0, ! ion positions and velocities
|
||||
> amass, ! atomic masses
|
||||
> katm, ! katm list, integer katm(nion)
|
||||
> katm_qm, ! katm list, integer katm(nion)
|
||||
> natm, ! natm list, integer natm(nkatm)
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@
|
|||
* ***************************
|
||||
subroutine xyz_init()
|
||||
*
|
||||
* $Id: xyz.F,v 1.8 2004-12-31 00:46:13 edo Exp $
|
||||
* $Id: xyz.F,v 1.9 2007-03-22 20:46:20 bylaska Exp $
|
||||
*
|
||||
implicit none
|
||||
|
||||
|
|
@ -93,10 +93,10 @@
|
|||
|
||||
* **** external functions ***
|
||||
character*2 ion_aname
|
||||
integer ion_nion, pspw_qmmm_nion
|
||||
integer ion_nion
|
||||
real*8 ion_rion, ion_vion
|
||||
external ion_aname
|
||||
external ion_nion, pspw_qmmm_nion
|
||||
external ion_nion
|
||||
external ion_rion, ion_vion
|
||||
integer ierr
|
||||
|
||||
|
|
@ -104,10 +104,10 @@
|
|||
|
||||
if (taskid.eq.MASTER) then
|
||||
#if defined(PSCALE) || defined(__crayx1)
|
||||
write(18,110,IOSTAT=ierr) (ion_nion()+pspw_qmmm_nion())
|
||||
write(18,110,IOSTAT=ierr) ion_nion()
|
||||
if(ierr.ne.0) backspace 18
|
||||
#else
|
||||
write(18,110) (ion_nion()+pspw_qmmm_nion())
|
||||
write(18,110) ion_nion()
|
||||
#endif
|
||||
do i=1,ion_nion()
|
||||
write(18,111) ion_aname(i),
|
||||
|
|
@ -118,7 +118,6 @@
|
|||
> ion_vion(2,i)*AACONV,
|
||||
> ion_vion(3,i)*AACONV
|
||||
end do
|
||||
call pspw_qmmm_PrintXYZ(18)
|
||||
call util_flush(18)
|
||||
end if
|
||||
110 format(I5/)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
/*
|
||||
$Id: qmmm_parse.c,v 1.2 2007-03-21 03:31:28 bylaska Exp $
|
||||
$Id: qmmm_parse.c,v 1.3 2007-03-22 20:46:20 bylaska Exp $
|
||||
*/
|
||||
|
||||
#include <math.h>
|
||||
|
|
@ -159,7 +159,7 @@ Integer *n3;
|
|||
/* define linear grid */
|
||||
nrl = 2001;
|
||||
rmax = 40.0;
|
||||
drl = rmax/((double) (nrl-1));
|
||||
drl = rmax/((double)(nrl-1));
|
||||
|
||||
fp = fopen(infile,"r+");
|
||||
w = get_word(fp);
|
||||
|
|
@ -225,7 +225,11 @@ Integer *n3;
|
|||
rr2 = rr1*rl[i];
|
||||
ttt = (rc1 - rr1);
|
||||
sss = (rc2 - rr2);
|
||||
pspl[0][i] = -Zion*(ttt/sss);
|
||||
/* l'Hopital */
|
||||
if (fabs(sss)<1.0e-9)
|
||||
pspl[0][i] = (-Zion/rc) * ((double) n_sigma)/((double) (n_sigma+1));
|
||||
else
|
||||
pspl[0][i] = -Zion*(ttt/sss);
|
||||
psil[0][i] = 0.0;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,12 @@
|
|||
*
|
||||
* $Id: paw_cpmd.F,v 1.22 2006-08-13 01:03:26 bylaska Exp $
|
||||
* $Id: paw_cpmd.F,v 1.23 2007-03-22 20:46:20 bylaska Exp $
|
||||
* $Log: not supported by cvs2svn $
|
||||
* Revision 1.22 2006/08/13 01:03:26 bylaska
|
||||
* Checking in code not include in 5.0 release.
|
||||
* A chain algorithm was added to Nose-Hoover thermostats.
|
||||
* Preliminary implementation of a processor group decomposition added to pspw, i.e. parallel decomposition is over fft grid and electrons.
|
||||
* ...EJB
|
||||
*
|
||||
* Revision 1.20 2006/02/11 02:50:46 bylaska
|
||||
* GGA's using 1st derivative formulas have been added in core part of PAW....EJB
|
||||
*
|
||||
|
|
@ -211,11 +217,9 @@
|
|||
real*8 control_rti,control_rte
|
||||
real*8 ion_vion,ion_com_ke,ion_ke
|
||||
real*8 ion_Temperature,ion_com_Temperature
|
||||
real*8 pspw_qmmm_Temperature
|
||||
external control_rti,control_rte
|
||||
external ion_vion,ion_com_ke,ion_ke
|
||||
external ion_Temperature,ion_com_Temperature
|
||||
external pspw_qmmm_Temperature
|
||||
|
||||
real*8 control_tole,control_tolc,control_tolr,ion_rion
|
||||
external control_tole,control_tolc,control_tolr,ion_rion
|
||||
|
|
@ -228,12 +232,10 @@
|
|||
integer control_it_in,control_it_out,control_gga,control_version
|
||||
integer control_ngrid,pack_nwave,control_lmax_multipole
|
||||
integer ion_nion,ion_natm,ion_katm,ion_nkatm
|
||||
integer pspw_qmmm_nion
|
||||
external pack_nwave_all
|
||||
external control_it_in,control_it_out,control_gga,control_version
|
||||
external control_ngrid,pack_nwave,control_lmax_multipole
|
||||
external ion_nion,ion_natm,ion_katm,ion_nkatm
|
||||
external pspw_qmmm_nion
|
||||
|
||||
character*12 control_boundry
|
||||
external control_boundry
|
||||
|
|
@ -245,8 +247,8 @@
|
|||
real*8 nwpw_timing,paw_psi_CheckOrtho
|
||||
external nwpw_timing,paw_psi_CheckOrtho
|
||||
|
||||
logical control_Nose,pspw_qmmm_found
|
||||
external control_Nose,pspw_qmmm_found
|
||||
logical control_Nose
|
||||
external control_Nose
|
||||
|
||||
integer Nose_Mchain,Nose_Nchain
|
||||
external Nose_Mchain,Nose_Nchain
|
||||
|
|
@ -559,7 +561,7 @@ c call init_paw_density(ispin,lmax)
|
|||
|
||||
|
||||
* **** initialize QM/MM ****
|
||||
call pspw_qmmm_init(rtdb)
|
||||
c call pspw_qmmm_init(rtdb)
|
||||
|
||||
* **** initialize dplot ****
|
||||
call dplot_iteration_init()
|
||||
|
|
@ -645,7 +647,6 @@ c call init_paw_density(ispin,lmax)
|
|||
write(6,1210) gx,gy,gz
|
||||
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
|
||||
write(6,1181)
|
||||
write(6,1192) (I,ion_aname(I),
|
||||
|
|
@ -653,7 +654,6 @@ c call init_paw_density(ispin,lmax)
|
|||
write(6,1200) vcx,vcy,vcz
|
||||
write(6,1210) vgx,vgy,vgz
|
||||
|
||||
call pspw_qmmm_Print2(6)
|
||||
|
||||
|
||||
|
||||
|
|
@ -1111,7 +1111,6 @@ c > dbl_mb(lmd(1)),ierr)
|
|||
write(6,1200) cx,cy,cz
|
||||
write(6,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
|
||||
write(6,1421)
|
||||
write(6,1192) (I,ion_aname(I),
|
||||
|
|
@ -1119,7 +1118,6 @@ c > dbl_mb(lmd(1)),ierr)
|
|||
write(6,1200) vcx,vcy,vcz
|
||||
write(6,1210) vgx,vgy,vgz
|
||||
|
||||
call pspw_qmmm_Print2(6)
|
||||
|
||||
|
||||
write(6,*)
|
||||
|
|
@ -1136,15 +1134,15 @@ c > dbl_mb(lmd(1)),ierr)
|
|||
write(6,1472) ion_ke(),ion_ke()/ion_nion()
|
||||
|
||||
|
||||
if (pspw_qmmm_found()) then
|
||||
write(6,1700)
|
||||
write(6,1701)
|
||||
write(6,1702) E(11)
|
||||
write(6,1703) E(12)
|
||||
write(6,1704) E(13)
|
||||
write(6,1705) E(14)
|
||||
write(6,1706) E(15)
|
||||
end if
|
||||
c if (pspw_qmmm_found()) then
|
||||
c write(6,1700)
|
||||
c write(6,1701)
|
||||
c write(6,1702) E(11)
|
||||
c write(6,1703) E(12)
|
||||
c write(6,1704) E(13)
|
||||
c write(6,1705) E(14)
|
||||
c write(6,1706) E(15)
|
||||
c end if
|
||||
|
||||
|
||||
if (control_Nose()) then
|
||||
|
|
@ -1159,12 +1157,6 @@ c > dbl_mb(lmd(1)),ierr)
|
|||
* **** write out Temperatures ****
|
||||
write(6,1491) eke
|
||||
eki0 = ion_Temperature()
|
||||
if (pspw_qmmm_found()) then
|
||||
eki1 =pspw_qmmm_Temperature()
|
||||
sum = ion_nion() + pspw_qmmm_nion() - 2.0d0
|
||||
eki0 = eki0*((ion_nion()-2.0d0)/sum)
|
||||
> + eki1*((pspw_qmmm_nion()-2.0d0)/sum)
|
||||
end if
|
||||
write(6,1480) eki0
|
||||
write(6,1490) ion_com_Temperature()
|
||||
|
||||
|
|
@ -1213,7 +1205,6 @@ c > dbl_mb(lmd(1)),ierr)
|
|||
|
||||
|
||||
!* **** write geometry to rtdb ****
|
||||
call pspw_qmmm_write(rtdb)
|
||||
call ion_write(rtdb)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,8 +1,11 @@
|
|||
*
|
||||
* $Id: paw_inner_loop_md.F,v 1.11 2006-02-11 02:50:46 bylaska Exp $
|
||||
* $Id: paw_inner_loop_md.F,v 1.12 2007-03-22 20:46:20 bylaska Exp $
|
||||
*
|
||||
|
||||
* $Log: not supported by cvs2svn $
|
||||
* Revision 1.11 2006/02/11 02:50:46 bylaska
|
||||
* GGA's using 1st derivative formulas have been added in core part of PAW....EJB
|
||||
*
|
||||
* Revision 1.10 2005/02/09 02:38:57 bylaska
|
||||
* ..............EJB
|
||||
*
|
||||
|
|
@ -158,24 +161,6 @@
|
|||
external Nose_e_energy,Nose_r_energy,Nose_sse,Nose_ssr
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found,pspw_charge_found
|
||||
real*8 pspw_qmmm_ke
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found,pspw_charge_found
|
||||
external pspw_qmmm_ke
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
logical dplot_iteration_check
|
||||
external dplot_iteration_check
|
||||
|
||||
|
|
@ -188,7 +173,7 @@
|
|||
deltamm = 0.0d0
|
||||
gga = control_gga()
|
||||
|
||||
field_exist = pspw_qmmm_found().or.pspw_charge_found()
|
||||
!field_exist = pspw_charge_found()
|
||||
c allow_translation = control_allow_translation()
|
||||
|
||||
|
||||
|
|
@ -292,7 +277,6 @@ c allow_translation = control_allow_translation()
|
|||
* *** skip ion_shift if newton step ***
|
||||
if (verlet) call ion_shift()
|
||||
if (nose) call Nose_Shift()
|
||||
if ((verlet).and.(pspw_qmmm_found())) call pspw_qmmm_shift()
|
||||
|
||||
|
||||
* ********************************
|
||||
|
|
|
|||
59
src/nwpw/paw/paw_basis/paw_small_basis.F
Normal file
59
src/nwpw/paw/paw_basis/paw_small_basis.F
Normal file
|
|
@ -0,0 +1,59 @@
|
|||
*
|
||||
* $Id: paw_small_basis.F,v 1.1 2007-03-22 20:46:20 bylaska Exp $
|
||||
*
|
||||
|
||||
!**************************************************
|
||||
!
|
||||
! Name: init_paw_small_basis
|
||||
!
|
||||
! Purpose
|
||||
!
|
||||
! Created: 7/11/2002
|
||||
!**************************************************
|
||||
subroutine init_paw_small_basis()
|
||||
implicit none
|
||||
|
||||
#include "paw_basis.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "errquit.fh"
|
||||
|
||||
* !*** local variables ***
|
||||
integer it, ir, nt,ngrid
|
||||
integer l_nb,l_ng
|
||||
integer i_nb,i_ng
|
||||
integer i_lo,i_hi,ii,mult_l
|
||||
integer max_mult_l
|
||||
integer max_ic
|
||||
logical ok
|
||||
double precision log_amesh,amesh
|
||||
character*255 fname
|
||||
|
||||
c !*** external functions ***
|
||||
integer ion_nkatm
|
||||
external ion_nkatm
|
||||
|
||||
* *** check for existence of _basis files ***
|
||||
nt = ion_nkatm()
|
||||
|
||||
* **** allocating memory from stack ****
|
||||
ok = MA_push_get(MT_INT,nt, "nb_tmp_array", l_nb,i_nb)
|
||||
ok = ok.and.
|
||||
> MA_push_get(MT_INT,nt, "ng_tmp_array", l_ng,i_ng)
|
||||
if (.not.ok)
|
||||
> call errquit('init_paw_small_basis:out of stack memory',0,MA_ERR)
|
||||
|
||||
do it = 1, nt
|
||||
ngrid = paw_basis_ngrid(it)
|
||||
write(*,*) "small basis, ngrid=",it,ngrid
|
||||
end do
|
||||
|
||||
|
||||
* **** deallocating memory from stack ****
|
||||
ok = MA_pop_stack(l_ng)
|
||||
ok = ok.and.MA_pop_stack(l_nb)
|
||||
if (.not.ok)
|
||||
> call errquit('init_paw_small_basis:error popping stack',1,MA_ERR)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: bfgsminimize.F,v 1.23 2006-10-13 01:43:58 bylaska Exp $
|
||||
c $Id: bfgsminimize.F,v 1.24 2007-03-22 20:46:20 bylaska Exp $
|
||||
c
|
||||
|
||||
* ************************************************************
|
||||
|
|
@ -45,8 +45,8 @@ c
|
|||
|
||||
logical value
|
||||
integer neall,npack1
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
!real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -83,22 +83,6 @@ c
|
|||
external linesearch
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found
|
||||
|
|
@ -299,36 +283,8 @@ c call dng_2to1()
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
if (pspw_qmmm_found()) then
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_e = psi_1v_field()
|
||||
e_qmmm_q = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
* **** get pspw_charge energies ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: bfgsminimize2.F,v 1.5 2006-08-13 01:03:26 bylaska Exp $
|
||||
c $Id: bfgsminimize2.F,v 1.6 2007-03-22 20:46:20 bylaska Exp $
|
||||
c
|
||||
|
||||
* ************************************************************
|
||||
|
|
@ -45,8 +45,8 @@ c
|
|||
|
||||
logical value
|
||||
integer neall,npack1
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
!real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 ehsic,exsic,phsic,pxsic,ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -86,22 +86,6 @@ c
|
|||
external linesearch
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found
|
||||
|
|
@ -262,36 +246,8 @@ c call dng_2to1()
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
if (pspw_qmmm_found()) then
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_e = psi_1v_field()
|
||||
e_qmmm_q = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
* **** get pspw_charge energies ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: bfgsminimize3.F,v 1.4 2006-08-13 01:03:26 bylaska Exp $
|
||||
c $Id: bfgsminimize3.F,v 1.5 2007-03-22 20:46:20 bylaska Exp $
|
||||
c
|
||||
|
||||
* ************************************************************
|
||||
|
|
@ -42,8 +42,8 @@ c
|
|||
|
||||
logical value
|
||||
integer neall,npack1
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
!real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -79,22 +79,6 @@ c
|
|||
external linesearch
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found
|
||||
|
|
@ -258,36 +242,8 @@ c call dng_2to1()
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
if (pspw_qmmm_found()) then
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_e = psi_1v_field()
|
||||
e_qmmm_q = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
* **** get pspw_charge energies ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: bybminimize.F,v 1.18 2006-04-08 00:25:48 bylaska Exp $
|
||||
c $Id: bybminimize.F,v 1.19 2007-03-22 20:46:20 bylaska Exp $
|
||||
c
|
||||
|
||||
* ************************************************************
|
||||
|
|
@ -46,8 +46,8 @@ c
|
|||
|
||||
logical value,precondition,done,stalled,deltav_bad(4),oprint
|
||||
integer n2ft3d
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
!real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -82,22 +82,6 @@ c
|
|||
external psi_1eorbit
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found,control_precondition,pspw_HFX
|
||||
|
|
@ -361,35 +345,8 @@ c
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
if (pspw_qmmm_found()) then
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_e = psi_1v_field()
|
||||
e_qmmm_q = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
* **** get pspw_charge energies ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: cgmdv5.F,v 1.3 2007-01-02 18:36:51 bylaska Exp $
|
||||
* $Id: cgmdv5.F,v 1.4 2007-03-22 20:46:20 bylaska Exp $
|
||||
*
|
||||
|
||||
***********************************************************************
|
||||
|
|
@ -105,11 +105,11 @@
|
|||
integer pack_nwave_all
|
||||
integer control_it_in,control_it_out,control_gga,control_version
|
||||
integer control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
|
||||
integer pspw_qmmm_nion,pspw_charge_nion
|
||||
integer pspw_charge_nion
|
||||
external pack_nwave_all
|
||||
external control_it_in,control_it_out,control_gga,control_version
|
||||
external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
|
||||
external pspw_qmmm_nion,pspw_charge_nion
|
||||
external pspw_charge_nion
|
||||
|
||||
character*12 control_boundry
|
||||
external control_boundry
|
||||
|
|
@ -121,9 +121,9 @@
|
|||
external psp_semicore
|
||||
external psp_rcore,psp_ncore,psp_rlocal
|
||||
|
||||
logical psi_initialize, psi_finalize, pspw_qmmm_found
|
||||
logical psi_initialize, psi_finalize
|
||||
integer psi_ispin,psi_ne,electron_count,linesearch_count
|
||||
external psi_initialize, psi_finalize, pspw_qmmm_found
|
||||
external psi_initialize, psi_finalize
|
||||
external psi_ispin,psi_ne,electron_count,linesearch_count
|
||||
real*8 cgsd_energy,cgsd_noit_energy
|
||||
external cgsd_energy,cgsd_noit_energy
|
||||
|
|
@ -141,8 +141,8 @@
|
|||
external control_ks_algorithm
|
||||
real*8 control_ks_alpha
|
||||
external control_ks_alpha
|
||||
logical control_print,control_qmmm,control_balance
|
||||
external control_print,control_qmmm,control_balance
|
||||
logical control_print,control_balance
|
||||
external control_print,control_balance
|
||||
integer control_mapping
|
||||
external control_mapping
|
||||
integer control_ks_maxit_orb,control_ks_maxit_orbs
|
||||
|
|
@ -210,7 +210,7 @@
|
|||
nwave = 0
|
||||
minimizer = control_minimizer()
|
||||
mapping = control_mapping()
|
||||
qmmm = control_qmmm()
|
||||
!qmmm = control_qmmm()
|
||||
|
||||
|
||||
* **** initialize D3dB data structure ****
|
||||
|
|
@ -346,7 +346,7 @@ c end if
|
|||
write(6,1121) control_boundry(),control_version()
|
||||
if (psi_ispin().eq.1) write(6,1130) "restricted"
|
||||
if (psi_ispin().eq.2) write(6,1130) "unrestricted"
|
||||
if (qmmm) write(6,1122)
|
||||
!if (qmmm) write(6,1122)
|
||||
IF (control_gga().eq.-1) THEN
|
||||
write(6,1131) 'off'
|
||||
ELSE IF (control_gga().eq.0) THEN
|
||||
|
|
@ -428,7 +428,6 @@ c end if
|
|||
write(6,1200) cx,cy,cz
|
||||
write(6,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
call pspw_charge_Print(6)
|
||||
|
||||
write(6,1181)
|
||||
|
|
@ -437,7 +436,6 @@ c end if
|
|||
write(6,1200) vcx,vcy,vcz
|
||||
write(6,1210) vgx,vgy,vgz
|
||||
|
||||
call pspw_qmmm_Print2(6)
|
||||
!end if
|
||||
|
||||
!write(6,1220) psi_ne(1),psi_ne(psi_ispin()),' (fourier space)'
|
||||
|
|
@ -576,7 +574,6 @@ c end if
|
|||
nion = ion_nion()
|
||||
if (pspw_charge_found().and.
|
||||
> (.not.pspw_bqext())) nion = nion + pspw_charge_nion()
|
||||
if (pspw_qmmm_found()) nion = nion + pspw_qmmm_nion()
|
||||
|
||||
value = MA_push_get(mt_dbl,(3*nion),
|
||||
> 'fion',fion(2),fion(1))
|
||||
|
|
@ -771,7 +768,6 @@ c call control_up_print()
|
|||
|
||||
* **** write geometry to rtdb ****
|
||||
call pspw_charge_write(rtdb)
|
||||
call pspw_qmmm_write(rtdb)
|
||||
call ion_write(rtdb)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: cgminimize.F,v 1.23 2007-03-21 03:31:28 bylaska Exp $
|
||||
c $Id: cgminimize.F,v 1.24 2007-03-22 20:46:21 bylaska Exp $
|
||||
c
|
||||
|
||||
* ************************************************************
|
||||
|
|
@ -45,8 +45,7 @@ c complex*16 G1(nfft3d,nemax)
|
|||
|
||||
logical value
|
||||
integer neall,npack1
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 e_lj,e_q
|
||||
real*8 ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -83,21 +82,11 @@ c complex*16 G1(nfft3d,nemax)
|
|||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
real*8 pspw_qmmm_LJ_E
|
||||
real*8 pspw_qmmm_Q_E
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_LJ_E
|
||||
external pspw_qmmm_Q_E
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found
|
||||
|
|
@ -316,30 +305,15 @@ c call dng_2to1()
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
|
||||
* **** set pspw_qmmm energies ****
|
||||
if (pspw_qmmm_found()) then
|
||||
e_lj = pspw_qmmm_LJ_E()
|
||||
e_q = pspw_qmmm_Q_E()
|
||||
E(1) = E(1) + e_lj + e_q
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
E(11) = e_lj
|
||||
E(12) = e_q
|
||||
end if
|
||||
|
||||
* **** get pspw_charge energies ****
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: cgminimize2.F,v 1.11 2006-10-07 00:10:05 bylaska Exp $
|
||||
c $Id: cgminimize2.F,v 1.12 2007-03-22 20:46:21 bylaska Exp $
|
||||
c
|
||||
|
||||
* **************************************************************
|
||||
|
|
@ -40,8 +40,8 @@ c
|
|||
|
||||
logical value
|
||||
integer neall,npack1
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
!real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 ehsic,exsic,phsic,pxsic,ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -77,22 +77,6 @@ c
|
|||
external linesearch
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found
|
||||
|
|
@ -264,36 +248,8 @@ c
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
if (pspw_qmmm_found()) then
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_e = psi_1v_field()
|
||||
e_qmmm_q = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
c $Id: cgsd_energy.F,v 1.66 2006-07-20 17:27:55 bylaska Exp $
|
||||
c $Id: cgsd_energy.F,v 1.67 2007-03-22 20:46:21 bylaska Exp $
|
||||
* ***************************
|
||||
* * *
|
||||
* * cgsd_energy *
|
||||
|
|
@ -36,7 +36,7 @@ c#include "rtdb.fh"
|
|||
integer rtdb
|
||||
|
||||
* **** external functions ****
|
||||
logical psp_semicore,pspw_qmmm_found,pspw_charge_found
|
||||
logical psp_semicore,pspw_charge_found,pspw_qmmm_found
|
||||
logical control_out_of_time,pspw_SIC,control_DOS,pspw_HFX
|
||||
logical control_print,control_psi_tmp
|
||||
character*2 ion_aname
|
||||
|
|
@ -45,7 +45,7 @@ c#include "rtdb.fh"
|
|||
integer psi_iptr_hml
|
||||
real*8 control_tole,control_tolc
|
||||
real*8 ion_rion,psi_eigenvalue,psi_hml,ion_amass
|
||||
external psp_semicore,pspw_qmmm_found,pspw_charge_found
|
||||
external psp_semicore,pspw_charge_found,pspw_qmmm_found
|
||||
external control_out_of_time,pspw_SIC,control_DOS,pspw_HFX
|
||||
external ion_aname
|
||||
external control_it_in, control_it_out,control_minimizer
|
||||
|
|
@ -229,11 +229,10 @@ c#include "rtdb.fh"
|
|||
if (pspw_qmmm_found()) then
|
||||
write(luout,1431)
|
||||
write(luout,1432)
|
||||
write(luout,1433)
|
||||
>(E(1)-E(11)-E(12)-E(13)-E(14)-E(15)-E(16)-E(17)-E(18)),
|
||||
>(E(1)-E(11)-E(12)-E(13)-E(14)-E(15)-E(16)-E(17)-E(18))
|
||||
> /ion_nion()
|
||||
write(luout,1433) E(1)-E(11),
|
||||
> (E(1)-E(11))/ion_nion()
|
||||
end if
|
||||
|
||||
if (pspw_charge_found()) then
|
||||
write(luout,1431)
|
||||
write(luout,1432)
|
||||
|
|
@ -254,11 +253,14 @@ c#include "rtdb.fh"
|
|||
write(luout,1470) E(5),E(5)/ion_nion()
|
||||
|
||||
write(luout,1480) E(6),E(6)/neall
|
||||
if (pspw_qmmm_found()) then
|
||||
write(luout,1491) E(7),E(7)/neall
|
||||
else
|
||||
write(luout,1490) E(7),E(7)/neall
|
||||
end if
|
||||
|
||||
c if (pspw_qmmm_found()) then
|
||||
c write(luout,1491) E(7),E(7)/neall
|
||||
c else
|
||||
c write(luout,1490) E(7),E(7)/neall
|
||||
c end if
|
||||
write(luout,1490) E(7),E(7)/neall
|
||||
|
||||
write(luout,1495) E(8),E(8)/neall
|
||||
write(luout,1496) E(9),E(9)/neall
|
||||
write(luout,1497) E(10),E(10)/neall
|
||||
|
|
@ -276,19 +278,17 @@ c#include "rtdb.fh"
|
|||
if (pspw_qmmm_found()) then
|
||||
write(luout,1700)
|
||||
write(luout,1701)
|
||||
write(luout,1702) E(11)+E(12)+E(13)
|
||||
write(luout,1707) E(11)
|
||||
write(luout,1708) E(12)
|
||||
write(luout,1709) E(13)
|
||||
|
||||
write(luout,1703) E(14)+E(15)
|
||||
write(luout,1710) E(14)
|
||||
write (luout,1711) E(15)
|
||||
|
||||
write(luout,1704) E(16)
|
||||
write(luout,1705) E(17)
|
||||
write(luout,1706) E(18)
|
||||
end if
|
||||
write(luout,1702) E(11)
|
||||
write(luout,1703) E(12)
|
||||
c
|
||||
c write(luout,1703) E(14)+E(15)
|
||||
c write(luout,1710) E(14)
|
||||
c write (luout,1711) E(15)
|
||||
c
|
||||
c write(luout,1704) E(16)
|
||||
c write(luout,1705) E(17)
|
||||
c write(luout,1706) E(18)
|
||||
end if
|
||||
if (pspw_charge_found()) then
|
||||
write(luout,1800)
|
||||
write(luout,1801)
|
||||
|
|
@ -472,8 +472,8 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
|
||||
1700 FORMAT(/' QM/MM-pol-vib/CAV Energies')
|
||||
1701 FORMAT( ' --------------------------')
|
||||
1702 FORMAT( ' QM/MM energy :',E19.10)
|
||||
1703 FORMAT(/' MM/MM energy :',E19.10)
|
||||
1702 FORMAT( ' LJ energy :',E19.10)
|
||||
1703 FORMAT(/' Residual Coulomb energy :',E19.10)
|
||||
1704 FORMAT(/' MM Polarization energy :',E19.10)
|
||||
1705 FORMAT( ' MM Vibration energy :',E19.10)
|
||||
1706 FORMAT( ' (QM+MM)/Cavity energy :',E19.10)
|
||||
|
|
@ -515,19 +515,19 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
real*8 fmx2,fmy2,fmz2
|
||||
|
||||
* **** external functions ****
|
||||
logical psp_semicore,pspw_qmmm_found,pspw_charge_found
|
||||
logical psp_semicore,pspw_charge_found,pspw_qmmm_found
|
||||
logical control_allow_translation,ion_q_FixIon,control_print
|
||||
character*2 ion_aname,pspw_qmmm_aname,pspw_charge_aname
|
||||
character*2 ion_aname,pspw_charge_aname
|
||||
integer ion_katm,ion_nion,control_version
|
||||
integer pspw_qmmm_nion,pspw_charge_nion
|
||||
real*8 ion_rion,pspw_qmmm_rion,pspw_charge_rion
|
||||
integer pspw_charge_nion
|
||||
real*8 ion_rion,pspw_charge_rion
|
||||
real*8 pspw_charge_charge
|
||||
external psp_semicore,pspw_qmmm_found,pspw_charge_found
|
||||
external psp_semicore,pspw_charge_found,pspw_qmmm_found
|
||||
external control_allow_translation,ion_q_FixIon,control_print
|
||||
external ion_aname,pspw_qmmm_aname,pspw_charge_aname
|
||||
external ion_aname,pspw_charge_aname
|
||||
external ion_katm,ion_nion,control_version
|
||||
external pspw_qmmm_nion,pspw_charge_nion
|
||||
external ion_rion,pspw_qmmm_rion,pspw_charge_rion
|
||||
external pspw_charge_nion
|
||||
external ion_rion,pspw_charge_rion
|
||||
external pspw_charge_charge
|
||||
logical pspw_bqext
|
||||
external pspw_bqext
|
||||
|
|
@ -537,7 +537,6 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
nion = ion_nion()
|
||||
if (pspw_charge_found().and.
|
||||
> (.not.pspw_bqext())) nion = nion + pspw_charge_nion()
|
||||
if (pspw_qmmm_found()) nion = nion + pspw_qmmm_nion()
|
||||
|
||||
call dcopy(3*nion,0.0d0,0,G1,1)
|
||||
|
||||
|
|
@ -546,6 +545,7 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
|
||||
if (control_version().eq.3) call ewald_f(G1)
|
||||
if (control_version().eq.4) call ion_ion_f(G1)
|
||||
if (pspw_qmmm_found()) call pspw_qmmm_fion(G1)
|
||||
|
||||
if (pspw_charge_found()) then
|
||||
if(pspw_bqext()) then
|
||||
|
|
@ -558,10 +558,6 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
end if
|
||||
end if
|
||||
|
||||
if (pspw_qmmm_found()) then
|
||||
nion1 = ion_nion() + pspw_charge_nion()
|
||||
call rho_1_Fion_Fmm(G1,G1(1,nion1+1))
|
||||
end if
|
||||
|
||||
* **** remove ion forces using ion_FixIon ****
|
||||
call ion_FixIon(G1)
|
||||
|
|
@ -625,21 +621,8 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
end if
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
|
||||
* **** print out mm positions ***
|
||||
if (pspw_qmmm_found()) then
|
||||
do i=1,pspw_qmmm_nion()
|
||||
i1 = ion_nion() + pspw_charge_nion() + i
|
||||
if (ion_q_FixIon(i1)) then
|
||||
write(luout,1191) i1,pspw_qmmm_aname(i),
|
||||
> (pspw_qmmm_rion(K,i),K=1,3)
|
||||
else
|
||||
write(luout,1190) i1,pspw_qmmm_aname(i),
|
||||
> (pspw_qmmm_rion(K,i),K=1,3)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end if
|
||||
|
||||
if (lprint) then
|
||||
write(luout,1421)
|
||||
|
|
@ -655,15 +638,6 @@ c1410 FORMAT(10X,'============= summary of results =================')
|
|||
end do
|
||||
end if
|
||||
|
||||
* **** print out mm forces ***
|
||||
if (pspw_qmmm_found()) then
|
||||
do i=1,pspw_qmmm_nion()
|
||||
i1 = ion_nion() + pspw_charge_nion() + i
|
||||
write(luout,1190) i1,pspw_qmmm_aname(i),
|
||||
> (G1(K,i1),K=1,3)
|
||||
end do
|
||||
end if
|
||||
|
||||
write(luout,1210) fmx2,fmy2,fmz2
|
||||
write(luout,1425)
|
||||
write(luout,1426) dsqrt(GG)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: cgsdv5.F,v 1.77 2007-03-19 19:05:36 bylaska Exp $
|
||||
* $Id: cgsdv5.F,v 1.78 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
|
||||
***********************************************************************
|
||||
|
|
@ -112,9 +112,9 @@
|
|||
external psp_semicore
|
||||
external psp_rcore,psp_ncore,psp_rlocal
|
||||
|
||||
logical psi_initialize, psi_finalize, pspw_qmmm_found
|
||||
logical psi_initialize, psi_finalize
|
||||
integer psi_ispin,psi_ne,psi_neq,electron_count,linesearch_count
|
||||
external psi_initialize, psi_finalize, pspw_qmmm_found
|
||||
external psi_initialize, psi_finalize
|
||||
external psi_ispin,psi_ne,psi_neq,electron_count,linesearch_count
|
||||
real*8 cgsd_energy,cgsd_noit_energy
|
||||
external cgsd_energy,cgsd_noit_energy
|
||||
|
|
@ -132,8 +132,8 @@
|
|||
external control_ks_algorithm
|
||||
real*8 control_ks_alpha
|
||||
external control_ks_alpha
|
||||
logical control_print,control_qmmm,control_balance
|
||||
external control_print,control_qmmm,control_balance
|
||||
logical control_print,control_balance
|
||||
external control_print,control_balance
|
||||
integer control_mapping,control_mapping1d,control_np_orbital
|
||||
external control_mapping,control_mapping1d,control_np_orbital
|
||||
integer control_ks_maxit_orb,control_ks_maxit_orbs
|
||||
|
|
@ -204,7 +204,6 @@
|
|||
nwave = 0
|
||||
minimizer = control_minimizer()
|
||||
mapping = control_mapping()
|
||||
qmmm = control_qmmm()
|
||||
|
||||
|
||||
* **** initialize D3dB data structure ****
|
||||
|
|
@ -315,7 +314,7 @@ c value = psi_initialize()
|
|||
write(luout,1121) control_boundry(),control_version()
|
||||
if (psi_ispin().eq.1) write(luout,1130) "restricted"
|
||||
if (psi_ispin().eq.2) write(luout,1130) "unrestricted"
|
||||
if (qmmm) write(luout,1122)
|
||||
!if (qmmm) write(luout,1122)
|
||||
IF (control_gga().eq.-1) THEN
|
||||
write(luout,1131) 'off'
|
||||
ELSE IF (control_gga().eq.0) THEN
|
||||
|
|
@ -397,7 +396,6 @@ c value = psi_initialize()
|
|||
write(luout,1200) cx,cy,cz
|
||||
write(luout,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(luout)
|
||||
call pspw_charge_Print(luout)
|
||||
end if
|
||||
|
||||
|
|
@ -580,7 +578,6 @@ c end if
|
|||
|
||||
* **** write geometry to rtdb ****
|
||||
call pspw_charge_write(rtdb)
|
||||
call pspw_qmmm_write(rtdb)
|
||||
call ion_write(rtdb)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: nominimize.F,v 1.6 2006-08-13 01:03:26 bylaska Exp $
|
||||
c $Id: nominimize.F,v 1.7 2007-03-22 20:46:21 bylaska Exp $
|
||||
c
|
||||
|
||||
subroutine nominimize(E)
|
||||
|
|
@ -19,8 +19,8 @@ c
|
|||
|
||||
logical value
|
||||
integer neall,npack1
|
||||
real*8 e_pol,e_vib,e_cav
|
||||
real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
!real*8 e_pol,e_vib,e_cav
|
||||
!real*8 e_qmmm_e,e_qmmm_q,e_qmmm_lj,e_mmmm_q,e_mmmm_lj
|
||||
real*8 ehsic,phsic,exsic,pxsic,ehfx,phfx
|
||||
|
||||
|
||||
|
|
@ -47,22 +47,6 @@ c
|
|||
external linesearch
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
* ***** pspw_charge external functions ****
|
||||
logical pspw_charge_found
|
||||
|
|
@ -110,36 +94,8 @@ c call psi_1toelectron()
|
|||
E(8) = psi_1vnl()
|
||||
E(9) = 2.0d0*ehartree
|
||||
E(10) = pxc
|
||||
if (pspw_qmmm_found()) then
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_qmmm_e = psi_1v_field()
|
||||
e_qmmm_q = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm_lj = pspw_qmmm_Energy_QMMM_LJ()
|
||||
|
||||
e_mmmm_q = pspw_qmmm_Energy_MMMM_Q()
|
||||
e_mmmm_lj = pspw_qmmm_Energy_MMMM_LJ()
|
||||
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_qmmm_q
|
||||
> + e_qmmm_lj
|
||||
> + e_mmmm_q
|
||||
> + e_mmmm_lj
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm_e
|
||||
E(12) = e_qmmm_q
|
||||
E(13) = e_qmmm_lj
|
||||
E(14) = e_mmmm_q
|
||||
E(15) = e_mmmm_lj
|
||||
E(16) = e_pol
|
||||
E(17) = e_vib
|
||||
E(18) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
* **** get pspw_charge energies ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: pspw_bq_gradient.F,v 1.1 2005-12-29 03:05:57 marat Exp $
|
||||
* $Id: pspw_bq_gradient.F,v 1.2 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
|
||||
|
||||
|
|
@ -88,11 +88,11 @@
|
|||
integer pack_nwave_all
|
||||
integer control_it_in,control_it_out,control_gga,control_version
|
||||
integer control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
|
||||
integer pspw_qmmm_nion,pspw_charge_nion
|
||||
integer pspw_charge_nion
|
||||
external pack_nwave_all
|
||||
external control_it_in,control_it_out,control_gga,control_version
|
||||
external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
|
||||
external pspw_qmmm_nion,pspw_charge_nion
|
||||
external pspw_charge_nion
|
||||
|
||||
character*12 control_boundry
|
||||
external control_boundry
|
||||
|
|
@ -104,9 +104,9 @@
|
|||
external psp_semicore
|
||||
external psp_rcore,psp_ncore,psp_rlocal
|
||||
|
||||
logical psi_initialize, psi_finalize, pspw_qmmm_found
|
||||
logical psi_initialize, psi_finalize
|
||||
integer psi_ispin,psi_ne,electron_count,linesearch_count
|
||||
external psi_initialize, psi_finalize, pspw_qmmm_found
|
||||
external psi_initialize, psi_finalize
|
||||
external psi_ispin,psi_ne,electron_count,linesearch_count
|
||||
real*8 cgsd_energy,cgsd_noit_energy
|
||||
external cgsd_energy,cgsd_noit_energy
|
||||
|
|
@ -124,8 +124,8 @@
|
|||
external control_ks_algorithm
|
||||
real*8 control_ks_alpha
|
||||
external control_ks_alpha
|
||||
logical control_print,control_qmmm,control_balance
|
||||
external control_print,control_qmmm,control_balance
|
||||
logical control_print,control_balance
|
||||
external control_print,control_balance
|
||||
integer control_mapping
|
||||
external control_mapping
|
||||
integer control_ks_maxit_orb,control_ks_maxit_orbs
|
||||
|
|
@ -189,7 +189,7 @@
|
|||
nwave = 0
|
||||
minimizer = control_minimizer()
|
||||
mapping = control_mapping()
|
||||
qmmm = control_qmmm()
|
||||
!qmmm = control_qmmm()
|
||||
|
||||
|
||||
* **** initialize D3dB data structure ****
|
||||
|
|
@ -269,8 +269,6 @@ c end if
|
|||
|
||||
* **** calculate gradient ***
|
||||
call dcopy(3*n_bq,0.0d0,0,f_bq,1)
|
||||
call pspw_charge_ion_Fcharge_qmmm(r_bq,q_bq,n_bq,f_bq)
|
||||
call rho_1Fcharge_qmmm(r_bq,q_bq,n_bq,f_bq)
|
||||
|
||||
if (taskid.eq.MASTER) call current_second(cpu3)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: cpmdv5.F,v 1.56 2007-02-23 01:24:32 bylaska Exp $
|
||||
* $Id: cpmdv5.F,v 1.57 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
|
||||
***********************************************************************
|
||||
|
|
@ -119,11 +119,9 @@ c parameter (autoatm =290.360032539d6)
|
|||
real*8 control_rti,control_rte,ion_rion
|
||||
real*8 ion_vion,ion_com_ke,ion_ke
|
||||
real*8 ion_Temperature,ion_com_Temperature
|
||||
real*8 pspw_qmmm_Temperature
|
||||
external control_rti,control_rte,ion_rion
|
||||
external ion_vion,ion_com_ke,ion_ke
|
||||
external ion_Temperature,ion_com_Temperature
|
||||
external pspw_qmmm_Temperature
|
||||
real*8 control_time_step,control_fake_mass
|
||||
external control_time_step,control_fake_mass
|
||||
logical control_read,control_move,ion_init,ion_q_FixIon
|
||||
|
|
@ -132,11 +130,9 @@ c parameter (autoatm =290.360032539d6)
|
|||
integer pack_nwave_all
|
||||
integer control_it_in,control_it_out,control_gga,control_version
|
||||
integer control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
|
||||
integer pspw_qmmm_nion
|
||||
external pack_nwave_all
|
||||
external control_it_in,control_it_out,control_gga,control_version
|
||||
external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm
|
||||
external pspw_qmmm_nion
|
||||
|
||||
character*12 control_boundry
|
||||
external control_boundry
|
||||
|
|
@ -146,8 +142,8 @@ c parameter (autoatm =290.360032539d6)
|
|||
external psp_semicore
|
||||
external psp_rcore,psp_ncore,psp_rlocal
|
||||
|
||||
logical control_Nose,control_Mulliken,pspw_qmmm_found
|
||||
external control_Nose,control_Mulliken,pspw_qmmm_found
|
||||
logical control_Nose,control_Mulliken
|
||||
external control_Nose,control_Mulliken
|
||||
|
||||
integer Nose_Mchain,Nose_Nchain
|
||||
external Nose_Mchain,Nose_Nchain
|
||||
|
|
@ -521,7 +517,6 @@ c parameter (autoatm =290.360032539d6)
|
|||
write(6,1200) cx,cy,cz
|
||||
write(6,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
|
||||
write(6,1181)
|
||||
write(6,1192) (I,ion_aname(I),
|
||||
|
|
@ -529,7 +524,6 @@ c parameter (autoatm =290.360032539d6)
|
|||
write(6,1200) vcx,vcy,vcz
|
||||
write(6,1210) vgx,vgy,vgz
|
||||
|
||||
call pspw_qmmm_Print2(6)
|
||||
|
||||
if (fractional) then
|
||||
write(6,1219) en(1),en(ispin),' ( fractional)'
|
||||
|
|
@ -989,7 +983,6 @@ c parameter (autoatm =290.360032539d6)
|
|||
write(6,1200) cx,cy,cz
|
||||
write(6,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
|
||||
write(6,1421)
|
||||
write(6,1192) (I,ion_aname(I),
|
||||
|
|
@ -997,7 +990,6 @@ c parameter (autoatm =290.360032539d6)
|
|||
write(6,1200) vcx,vcy,vcz
|
||||
write(6,1210) vgx,vgy,vgz
|
||||
|
||||
call pspw_qmmm_Print2(6)
|
||||
|
||||
write(6,*)
|
||||
write(6,1320) en(1),en(ispin),' (real space)'
|
||||
|
|
@ -1017,15 +1009,15 @@ c parameter (autoatm =290.360032539d6)
|
|||
write(6,1472) ion_ke(),ion_ke()/ion_nion()
|
||||
|
||||
|
||||
if (pspw_qmmm_found()) then
|
||||
write(6,1700)
|
||||
write(6,1701)
|
||||
write(6,1702) E(11)
|
||||
write(6,1703) E(12)
|
||||
write(6,1704) E(13)
|
||||
write(6,1705) E(14)
|
||||
write(6,1706) E(15)
|
||||
end if
|
||||
c if (pspw_qmmm_found()) then
|
||||
c write(6,1700)
|
||||
c write(6,1701)
|
||||
c write(6,1702) E(11)
|
||||
c write(6,1703) E(12)
|
||||
c write(6,1704) E(13)
|
||||
c write(6,1705) E(14)
|
||||
c write(6,1706) E(15)
|
||||
c end if
|
||||
|
||||
|
||||
if (control_Nose()) then
|
||||
|
|
@ -1040,12 +1032,12 @@ c parameter (autoatm =290.360032539d6)
|
|||
* **** write out Temperatures ****
|
||||
write(6,1491) eke
|
||||
eki0 = ion_Temperature()
|
||||
if (pspw_qmmm_found()) then
|
||||
eki1 =pspw_qmmm_Temperature()
|
||||
sum = ion_nion() + pspw_qmmm_nion() - 2.0d0
|
||||
eki0 = eki0*((ion_nion()-2.0d0)/sum)
|
||||
> + eki1*((pspw_qmmm_nion()-2.0d0)/sum)
|
||||
end if
|
||||
c if (pspw_qmmm_found()) then
|
||||
c eki1 =pspw_qmmm_Temperature()
|
||||
c sum = ion_nion() + pspw_qmmm_nion() - 2.0d0
|
||||
c eki0 = eki0*((ion_nion()-2.0d0)/sum)
|
||||
c > + eki1*((pspw_qmmm_nion()-2.0d0)/sum)
|
||||
c end if
|
||||
write(6,1480) eki0
|
||||
write(6,1490) ion_com_Temperature()
|
||||
|
||||
|
|
@ -1128,7 +1120,6 @@ c parameter (autoatm =290.360032539d6)
|
|||
call v_psi_write(ispin,ne,dcpl_mb(psi0(1)))
|
||||
|
||||
* **** write geometry to rtdb ****
|
||||
call pspw_qmmm_write(rtdb)
|
||||
call ion_write(rtdb)
|
||||
|
||||
* **** deallocate heap memory ****
|
||||
|
|
|
|||
|
|
@ -14,7 +14,7 @@
|
|||
|
||||
logical function cpsdv5(rtdb)
|
||||
*
|
||||
* $Id: cpsdv5.F,v 1.49 2007-03-19 19:05:37 bylaska Exp $
|
||||
* $Id: cpsdv5.F,v 1.50 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
implicit none
|
||||
integer rtdb
|
||||
|
|
@ -131,13 +131,13 @@
|
|||
character*12 control_boundry
|
||||
external control_boundry
|
||||
|
||||
logical pspw_SIC,pspw_SIC_relaxed
|
||||
logical pspw_SIC,pspw_SIC_relaxed,pspw_qmmm_found
|
||||
logical pspw_HFX,pspw_HFX_relaxed
|
||||
logical psp_semicore,control_Mulliken,pspw_qmmm_found
|
||||
logical psp_semicore,control_Mulliken
|
||||
real*8 psp_rcore,psp_ncore,psp_rlocal
|
||||
external pspw_SIC,pspw_SIC_relaxed
|
||||
external pspw_SIC,pspw_SIC_relaxed,pspw_qmmm_found
|
||||
external pspw_HFX,pspw_HFX_relaxed
|
||||
external psp_semicore,control_Mulliken,pspw_qmmm_found
|
||||
external psp_semicore,control_Mulliken
|
||||
external psp_rcore,psp_ncore,psp_rlocal
|
||||
logical control_check_charge_multiplicity
|
||||
external control_check_charge_multiplicity
|
||||
|
|
@ -491,7 +491,6 @@ c end if
|
|||
write(6,1200) cx,cy,cz
|
||||
write(6,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
|
||||
if (fractional) then
|
||||
write(6,1219) en(1),en(ispin),' ( fractional)'
|
||||
|
|
@ -700,7 +699,6 @@ c end if
|
|||
write(6,1200) cx,cy,cz
|
||||
write(6,1210) gx,gy,gz
|
||||
|
||||
call pspw_qmmm_Print(6)
|
||||
|
||||
write(6,*)
|
||||
write(6,1320) en(1),en(ispin),' (real space)'
|
||||
|
|
@ -709,8 +707,8 @@ c end if
|
|||
if (pspw_qmmm_found()) then
|
||||
write(6,1431)
|
||||
write(6,1432)
|
||||
write(6,1433) (E(1)-E(11)-E(12)-E(13)-E(14)-E(15)),
|
||||
> (E(1)-E(11)-E(12)-E(13)-E(14)-E(15))/ion_nion()
|
||||
write(6,1433) (E(1)-E(11)),
|
||||
> (E(1)-E(11))/ion_nion()
|
||||
end if
|
||||
write(6,1440) E(2),E(2)/n2(ispin)
|
||||
write(6,1450) E(3),E(3)/n2(ispin)
|
||||
|
|
@ -739,15 +737,21 @@ c end if
|
|||
virial = (E(10)+E(9)+E(8)+E(7))/E(6)
|
||||
write(6,1498) virial
|
||||
|
||||
if (pspw_qmmm_found()) then
|
||||
write(6,1700)
|
||||
write(6,1701)
|
||||
write(6,1702) E(11)
|
||||
if (pspw_qmmm_found()) then
|
||||
write(6,1700)
|
||||
write(6,1701)
|
||||
write(6,1702) E(11)
|
||||
write(6,1703) E(12)
|
||||
write(6,1704) E(13)
|
||||
write(6,1705) E(14)
|
||||
write(6,1706) E(15)
|
||||
end if
|
||||
c
|
||||
c write(luout,1703) E(14)+E(15)
|
||||
c write(luout,1710) E(14)
|
||||
c write (luout,1711) E(15)
|
||||
c
|
||||
c write(luout,1704) E(16)
|
||||
c write(luout,1705) E(17)
|
||||
c write(luout,1706) E(18)
|
||||
end if
|
||||
|
||||
|
||||
* **** write out diagonal <psi|H|psi> matrix ****
|
||||
if (pspw_SIC()) then
|
||||
|
|
@ -832,7 +836,6 @@ c end if
|
|||
call psi_write(ispin,ne,dcpl_mb(psi2(1)))
|
||||
|
||||
* **** write geometry to rtdb ****
|
||||
call pspw_qmmm_write(rtdb)
|
||||
call ion_write(rtdb)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*,eke
|
||||
* $Id: inner_loop.F,v 1.24 2007-02-23 01:24:32 bylaska Exp $
|
||||
* $Id: inner_loop.F,v 1.25 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
|
||||
* ************************************************************
|
||||
|
|
@ -53,7 +53,8 @@
|
|||
real*8 exc,exc2,pxc,pxc2,dte,scal1,scal2,dv,dt
|
||||
real*8 ehsic,phsic,exsic,pxsic,ehfx,phfx
|
||||
real*8 deltamm
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 e_lj, e_q
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
|
||||
|
||||
|
||||
|
|
@ -102,29 +103,16 @@ c real*8 fion(3,natmx)
|
|||
external Dneall_m_push_get_block,Dneall_m_pop_stack
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found,pspw_charge_found
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found,pspw_charge_found
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
|
||||
logical pspw_charge_found,pspw_qmmm_found
|
||||
real*8 pspw_qmmm_LJ_E,pspw_qmmm_Q_E
|
||||
external pspw_charge_found,pspw_qmmm_found
|
||||
external pspw_qmmm_LJ_E,pspw_qmmm_Q_E
|
||||
|
||||
call Pack_npack(0,npack0)
|
||||
n2ft3d = 2*nfft3d
|
||||
deltamm = 0.0d0
|
||||
|
||||
field_exist = pspw_charge_found().or.pspw_qmmm_found()
|
||||
field_exist = pspw_charge_found()
|
||||
sic = pspw_SIC()
|
||||
hfx = pspw_HFX()
|
||||
gga = control_gga()
|
||||
|
|
@ -210,7 +198,6 @@ c real*8 fion(3,natmx)
|
|||
do it=1,it_in
|
||||
call dcopy(2*npack1*nemax,psi2,1,psi1,1)
|
||||
if (move) call ion_shift()
|
||||
if ((move).and.(pspw_qmmm_found())) call pspw_qmmm_shift()
|
||||
|
||||
* ********************************
|
||||
* **** generate phaze factors ****
|
||||
|
|
@ -356,12 +343,6 @@ c call mask_C(0,dcpl_mb(dng(1)))
|
|||
|
||||
call dcopy(n2ft3d,0.0d0,0,dbl_mb(v_field(1)),1)
|
||||
|
||||
* **** generate QM/MM potential ****
|
||||
if (pspw_qmmm_found()) then
|
||||
call pspw_qmmm_Generate_V(n2ft3d,
|
||||
> dbl_mb(r_grid(1)),
|
||||
> dbl_mb(v_field(1)))
|
||||
end if
|
||||
|
||||
* **** generate charge potential ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
@ -419,14 +400,8 @@ c call mask_C(0,dcpl_mb(dng(1)))
|
|||
call semicore_xc_F(ispin,dbl_mb(xcp(1)),dbl_mb(fion(1)))
|
||||
end if
|
||||
|
||||
* **** QM/MM update ****
|
||||
if (pspw_qmmm_found()) then
|
||||
call pspw_qmmm_Update(0,n2ft3d,dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),
|
||||
> dv,dt,
|
||||
> dbl_mb(fion(1)),
|
||||
> deltamm)
|
||||
end if
|
||||
* **** get the qmmm force ****
|
||||
if (pspw_qmmm_found()) call pspw_qmmm_fion(dbl_mb(fion(1)))
|
||||
|
||||
* **** remove ion forces using ion_FixIon ****
|
||||
call ion_FixIon(dbl_mb(fion(1)))
|
||||
|
|
@ -560,27 +535,12 @@ c > lmd,dbl_mb(tmp_L(1)),ierr)
|
|||
E(10) = pxc
|
||||
|
||||
if (pspw_qmmm_found()) then
|
||||
e_lj = pspw_qmmm_LJ_E()
|
||||
e_q = pspw_qmmm_Q_E()
|
||||
E(1) = E(1) + e_lj + e_q
|
||||
|
||||
* **** get QM/MM energies ****
|
||||
e_ionmm = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm = pspw_qmmm_Energy_QMMM_LJ()
|
||||
e_mmmm = pspw_qmmm_Energy_MMMM_Q()
|
||||
> + pspw_qmmm_Energy_MMMM_LJ()
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(1) = E(1) + e_ionmm
|
||||
> + e_qmmm
|
||||
> + e_mmmm
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(11) = e_qmmm
|
||||
E(12) = e_mmmm
|
||||
E(13) = e_pol
|
||||
E(14) = e_vib
|
||||
E(15) = e_cav
|
||||
E(11) = e_lj
|
||||
E(12) = e_q
|
||||
end if
|
||||
|
||||
* **** SIC corrections ****
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: inner_loop_md.F,v 1.26 2007-02-23 01:24:32 bylaska Exp $
|
||||
* $Id: inner_loop_md.F,v 1.27 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
* ************************************************************
|
||||
* * MPI cpmd routine *
|
||||
|
|
@ -58,7 +58,7 @@
|
|||
real*8 sum,Eold,eorbit,eion,ehartr,eke,eki,sse,ssr,sa1,sa2
|
||||
real*8 exc,exc2,pxc,pxc2,dte,dte0,scal1,scal2,dv,dt,fmass,h
|
||||
real*8 ehsic,phsic,exsic,pxsic,ehfx,phfx,espring,enlocal
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
!real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
real*8 s,r
|
||||
|
||||
|
||||
|
|
@ -103,24 +103,9 @@ c real*8 fion(3,natmx)
|
|||
external Nose_dXe,Nose_dXr
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
logical pspw_qmmm_found,pspw_charge_found
|
||||
real*8 pspw_qmmm_ke
|
||||
real*8 pspw_qmmm_Energy_QMMM_Q
|
||||
real*8 pspw_qmmm_Energy_QMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_MMMM_Q
|
||||
real*8 pspw_qmmm_Energy_MMMM_LJ
|
||||
real*8 pspw_qmmm_Energy_pol
|
||||
real*8 pspw_qmmm_Energy_vib
|
||||
real*8 pspw_qmmm_Energy_cav
|
||||
external pspw_qmmm_found,pspw_charge_found
|
||||
external pspw_qmmm_ke
|
||||
external pspw_qmmm_Energy_QMMM_Q
|
||||
external pspw_qmmm_Energy_QMMM_LJ
|
||||
external pspw_qmmm_Energy_MMMM_Q
|
||||
external pspw_qmmm_Energy_MMMM_LJ
|
||||
external pspw_qmmm_Energy_pol
|
||||
external pspw_qmmm_Energy_vib
|
||||
external pspw_qmmm_Energy_cav
|
||||
logical pspw_charge_found
|
||||
external pspw_charge_found
|
||||
|
||||
logical dplot_iteration_check
|
||||
external dplot_iteration_check
|
||||
|
||||
|
|
@ -138,7 +123,7 @@ c real*8 fion(3,natmx)
|
|||
call Pack_npack(0,npack0)
|
||||
|
||||
n2ft3d = 2*nfft3d
|
||||
field_exist = pspw_qmmm_found().or.pspw_charge_found()
|
||||
field_exist = pspw_charge_found()
|
||||
sic = pspw_SIC()
|
||||
gga = control_gga()
|
||||
allow_translation = control_allow_translation()
|
||||
|
|
@ -239,7 +224,7 @@ c real*8 fion(3,natmx)
|
|||
* *** skip ion_shift if newton step ***
|
||||
if (verlet) call ion_shift()
|
||||
if (nose.and.verlet) call Nose_Shift()
|
||||
if ((verlet).and.(pspw_qmmm_found())) call pspw_qmmm_shift()
|
||||
|
||||
|
||||
* ********************************
|
||||
* **** generate phaze factors ****
|
||||
|
|
@ -382,12 +367,6 @@ c real*8 fion(3,natmx)
|
|||
|
||||
call dcopy(n2ft3d,0.0d0,0,dbl_mb(v_field(1)),1)
|
||||
|
||||
* **** generate QM/MM potential ****
|
||||
if (pspw_qmmm_found()) then
|
||||
call pspw_qmmm_Generate_V(n2ft3d,
|
||||
> dbl_mb(r_grid(1)),
|
||||
> dbl_mb(v_field(1)))
|
||||
end if
|
||||
|
||||
* **** generate charge potential ****
|
||||
if (pspw_charge_found()) then
|
||||
|
|
@ -469,13 +448,6 @@ c end if
|
|||
end do
|
||||
|
||||
* **** QM/MM Verlet update ****
|
||||
if (pspw_qmmm_found()) then
|
||||
call pspw_qmmm_Update(2,n2ft3d,dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),
|
||||
> dv,dt,
|
||||
> dbl_mb(fion(1)),
|
||||
> sum)
|
||||
end if
|
||||
call ion_verlet_step(dbl_mb(fion(1)),sa_alpha(2))
|
||||
end if
|
||||
|
||||
|
|
@ -496,13 +468,6 @@ c end if
|
|||
end do
|
||||
|
||||
* **** QM/MM Newton update ****
|
||||
if (pspw_qmmm_found()) then
|
||||
call pspw_qmmm_Update(1,n2ft3d,dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),
|
||||
> dv,dt,
|
||||
> dbl_mb(fion(1)),
|
||||
> sum)
|
||||
end if
|
||||
call ion_newton_step(dbl_mb(fion(1)),sa_alpha(2)*r)
|
||||
|
||||
end if
|
||||
|
|
@ -667,28 +632,8 @@ c end if
|
|||
E(22) = espring
|
||||
|
||||
* ******** QM/MM energies ******
|
||||
if (pspw_qmmm_found()) then
|
||||
e_ionmm = pspw_qmmm_Energy_QMMM_Q()
|
||||
e_qmmm = pspw_qmmm_Energy_QMMM_LJ()
|
||||
e_mmmm = pspw_qmmm_Energy_MMMM_Q()
|
||||
> + pspw_qmmm_Energy_MMMM_LJ()
|
||||
e_pol = pspw_qmmm_Energy_pol()
|
||||
e_vib = pspw_qmmm_Energy_vib()
|
||||
e_cav = pspw_qmmm_Energy_cav()
|
||||
|
||||
E(2) = E(2) + e_ionmm
|
||||
> + e_qmmm
|
||||
> + e_mmmm
|
||||
> + e_pol
|
||||
> + e_vib
|
||||
> + e_cav
|
||||
E(4) = E(4) + pspw_qmmm_ke()
|
||||
E(11) = e_qmmm
|
||||
E(12) = e_mmmm
|
||||
E(13) = e_pol
|
||||
E(14) = e_vib
|
||||
E(15) = e_cav
|
||||
end if
|
||||
!if (pspw_qmmm_found()) then
|
||||
!end if
|
||||
|
||||
* **** SIC corrections ****
|
||||
if (pspw_SIC()) then
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: inner_loop_qmmm_step.F,v 1.1 2005-08-12 23:08:43 bylaska Exp $
|
||||
* $Id: inner_loop_qmmm_step.F,v 1.2 2007-03-22 20:46:21 bylaska Exp $
|
||||
*
|
||||
|
||||
subroutine inner_loop_qmmm_step(verlet,
|
||||
|
|
@ -47,7 +47,7 @@
|
|||
real*8 sum,Eold,eorbit,eion,ehartr,eke,eki,sse,ssr
|
||||
real*8 exc,exc2,pxc,pxc2,dte,dte0,scal1,scal2,dv,dt,fmass,h
|
||||
real*8 ehsic,phsic,exsic,pxsic,ehfx,phfx
|
||||
real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
c real*8 e_ionmm,e_qmmm,e_mmmm,e_pol,e_vib,e_cav
|
||||
|
||||
|
||||
* **** MA local variables ****
|
||||
|
|
@ -79,8 +79,8 @@
|
|||
external Nose_e_energy,Nose_r_energy,Nose_sse,Nose_ssr
|
||||
|
||||
* ***** QM/MM external functions ****
|
||||
real*8 pspw_charge_Energy_ion_qmmm
|
||||
external pspw_charge_Energy_ion_qmmm
|
||||
c real*8 pspw_charge_Energy_ion_qmmm
|
||||
c external pspw_charge_Energy_ion_qmmm
|
||||
logical dplot_iteration_check
|
||||
external dplot_iteration_check
|
||||
|
||||
|
|
@ -324,11 +324,11 @@
|
|||
|
||||
|
||||
* **** generate charge potential ****
|
||||
call pspw_charge_Generate_V_qmmm(rcharge,qcharge,ncharge,
|
||||
> n2ft3d,
|
||||
> dbl_mb(r_grid(1)),
|
||||
> dbl_mb(v_field(1)))
|
||||
|
||||
c call pspw_charge_Generate_V_qmmm(rcharge,qcharge,ncharge,
|
||||
c > n2ft3d,
|
||||
c > dbl_mb(r_grid(1)),
|
||||
c > dbl_mb(v_field(1)))
|
||||
c
|
||||
end if
|
||||
|
||||
|
||||
|
|
@ -370,14 +370,14 @@
|
|||
* **** get qmmm forces ****
|
||||
* **************************
|
||||
|
||||
call pspw_charge_rho_Fcharge_qmmm(rcharge,qcharge,ncharge,
|
||||
> n2ft3d,
|
||||
> dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),
|
||||
> dv,
|
||||
> fcharge)
|
||||
call pspw_charge_Fion_Fcharge_qmmm(rcharge,qcharge,ncharge,
|
||||
> fion,fcharge)
|
||||
c call pspw_charge_rho_Fcharge_qmmm(rcharge,qcharge,ncharge,
|
||||
c > n2ft3d,
|
||||
c > dbl_mb(r_grid(1)),
|
||||
c > dbl_mb(rho(1)),
|
||||
c > dv,
|
||||
c > fcharge)
|
||||
c call pspw_charge_Fion_Fcharge_qmmm(rcharge,qcharge,ncharge,
|
||||
c > fion,fcharge)
|
||||
|
||||
|
||||
* *****************************************
|
||||
|
|
@ -587,12 +587,12 @@
|
|||
if (ispin.eq.1) eke = 2.0d0*eke
|
||||
|
||||
|
||||
e_ionmm = pspw_charge_Energy_ion_qmmm(rcharge,qcharge,ncharge)
|
||||
c e_ionmm = pspw_charge_Energy_ion_qmmm(rcharge,qcharge,ncharge)
|
||||
|
||||
|
||||
* **** total energy ****
|
||||
Eold=E(1)
|
||||
E(2) = eorbit + eion + exc - ehartr - pxc + e_ionmm
|
||||
E(2) = eorbit + eion + exc - ehartr - pxc
|
||||
E(3) = eke
|
||||
E(4) = 0.0d0
|
||||
E(5) = eorbit
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: electron.F,v 1.48 2007-02-23 01:24:32 bylaska Exp $
|
||||
* $Id: electron.F,v 1.49 2007-03-22 20:46:22 bylaska Exp $
|
||||
*
|
||||
|
||||
|
||||
|
|
@ -1952,9 +1952,9 @@ c call Grsm_gg_Copy(npack1,nemax,psi_k,dbl_mb(psi_r(1)))
|
|||
integer r_grid(2)
|
||||
|
||||
* **** external functions *****
|
||||
logical pspw_qmmm_found,pspw_charge_found
|
||||
logical pspw_charge_found
|
||||
integer control_version
|
||||
external pspw_qmmm_found,pspw_charge_found
|
||||
external pspw_charge_found
|
||||
external control_version
|
||||
|
||||
value = MA_push_get(mt_dcpl,(nfft3d),'tmp1',tmp1(2),tmp1(1))
|
||||
|
|
@ -1980,7 +1980,6 @@ c call Grsm_gg_Copy(npack1,nemax,psi_k,dbl_mb(psi_r(1)))
|
|||
|
||||
* **** generate real-space fields ****
|
||||
if ((control_version().eq.4).or.
|
||||
> (pspw_qmmm_found()).or.
|
||||
> (pspw_charge_found())) then
|
||||
|
||||
value = MA_push_get(mt_dbl,(6*nfft3d),'r_grid',
|
||||
|
|
@ -1998,14 +1997,6 @@ c call Grsm_gg_Copy(npack1,nemax,psi_k,dbl_mb(psi_r(1)))
|
|||
call dcopy(2*nfft3d,0.0d0,0,dbl_mb(v_field(1)),1)
|
||||
|
||||
|
||||
* **** generate QM/MM potential ****
|
||||
if (pspw_qmmm_found()) then
|
||||
field_exist = .true.
|
||||
n2ft3d = 2*nfft3d
|
||||
call pspw_qmmm_Generate_V(n2ft3d,dbl_mb(r_grid(1)),
|
||||
> dbl_mb(v_field(1)))
|
||||
end if
|
||||
|
||||
* **** generate charge potential ****
|
||||
if (pspw_charge_found()) then
|
||||
field_exist = .true.
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: ewald.F,v 1.8 2006-10-07 00:10:06 bylaska Exp $
|
||||
* $Id: ewald.F,v 1.9 2007-03-22 20:46:22 bylaska Exp $
|
||||
*
|
||||
integer function ewald_ncut()
|
||||
implicit none
|
||||
|
|
@ -42,6 +42,19 @@
|
|||
return
|
||||
end
|
||||
|
||||
integer function ewald_rcell_ptr()
|
||||
implicit none
|
||||
|
||||
* **** common block for ewald.f ****
|
||||
integer ncut
|
||||
real*8 rcut,cewald,alpha
|
||||
integer vg(2),rcell(2)
|
||||
common / ewald_block / vg,rcell,cewald,alpha,rcut,ncut
|
||||
|
||||
ewald_rcell_ptr = rcell(1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
real*8 function ewald_mandelung()
|
||||
implicit none
|
||||
|
|
@ -369,7 +382,6 @@ c > + i+1
|
|||
real*8 function ewald_e()
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,9 +1,12 @@
|
|||
*
|
||||
* $Id: psi.F,v 1.75 2007-02-23 01:24:32 bylaska Exp $
|
||||
* $Id: psi.F,v 1.76 2007-03-22 20:46:22 bylaska Exp $
|
||||
*
|
||||
|
||||
|
||||
* $Log: not supported by cvs2svn $
|
||||
* Revision 1.75 2007/02/23 01:24:32 bylaska
|
||||
* ...EJB
|
||||
*
|
||||
* Revision 1.74 2007/02/10 03:56:54 bylaska
|
||||
* ...bug fix...
|
||||
* ..EJB
|
||||
|
|
@ -2624,73 +2627,6 @@ c call OrthoCheck(ispin,ne,dcpl_mb(psi1(1)))
|
|||
end
|
||||
|
||||
|
||||
* ***********************************
|
||||
* * *
|
||||
* * rho_1_Fion_Fmm *
|
||||
* * *
|
||||
* ***********************************
|
||||
subroutine rho_1_Fion_Fmm(fion,fmm)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
real*8 fion(*)
|
||||
real*8 fmm(*)
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "psi.fh"
|
||||
|
||||
* **** local variables ****
|
||||
logical value
|
||||
integer n2ft3d,nx,ny,nz
|
||||
integer r_grid(2),rho(2)
|
||||
real*8 dv
|
||||
|
||||
* **** external functions ****
|
||||
real*8 lattice_omega
|
||||
external lattice_omega
|
||||
|
||||
* **** Initializationsr ****
|
||||
call D3dB_n2ft3d(1,n2ft3d)
|
||||
call D3dB_nx(1,nx)
|
||||
call D3dB_ny(1,ny)
|
||||
call D3dB_nz(1,nz)
|
||||
dv = lattice_omega()/dble(nx*ny*nz)
|
||||
|
||||
|
||||
* **** Push memory ****
|
||||
value = MA_push_get(mt_dbl,(3*n2ft3d),'r_grid',
|
||||
> r_grid(2),r_grid(1))
|
||||
value = value.and.
|
||||
> MA_push_get(mt_dbl,(3*n2ft3d),'rho',
|
||||
> rho(2),rho(1))
|
||||
if (.not. value) call errquit(
|
||||
> 'rho_1_Fion_Fmm: out of stack memory',0, MA_ERR)
|
||||
|
||||
|
||||
* **** Get r_grid and rho ****
|
||||
call lattice_r_grid(dbl_mb(r_grid(1)))
|
||||
call D3dB_rr_Sum(1,dbl_mb(rho1(1)),
|
||||
> dbl_mb(rho1(1)+(ispin-1)*n2ft3d),
|
||||
> dbl_mb(rho(1)))
|
||||
|
||||
* **** Now calculate fion ****
|
||||
call pspw_qmmm_Fion(n2ft3d,dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),
|
||||
> dv,fion)
|
||||
|
||||
* **** Now calculate fmm ****
|
||||
call pspw_qmmm_Fmm(n2ft3d,dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),
|
||||
> dv,fmm)
|
||||
|
||||
* **** Pop memory ****
|
||||
value = MA_pop_stack(rho(2))
|
||||
value = value.and.MA_pop_stack(r_grid(2))
|
||||
if (.not. value) call errquit(
|
||||
> 'rho_1_Fion_Fmm: error popping stack memory',0, MA_ERR)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
* ***********************************
|
||||
* * *
|
||||
* * rho_1Fcharge *
|
||||
|
|
@ -2753,76 +2689,6 @@ c call OrthoCheck(ispin,ne,dcpl_mb(psi1(1)))
|
|||
end
|
||||
|
||||
|
||||
* ***********************************
|
||||
* * *
|
||||
* * rho_1Fcharge_qmmm *
|
||||
* * *
|
||||
* ***********************************
|
||||
subroutine rho_1Fcharge_qmmm(r_charge,q_charge,n_charge,f_charge)
|
||||
implicit none
|
||||
double precision r_charge(*)
|
||||
double precision q_charge(*)
|
||||
integer n_charge
|
||||
double precision f_charge(*)
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "psi.fh"
|
||||
#include "errquit.fh"
|
||||
|
||||
* **** local variables ****
|
||||
logical value
|
||||
integer n2ft3d,nx,ny,nz
|
||||
integer r_grid(2),rho(2)
|
||||
real*8 dv
|
||||
|
||||
* **** external functions ****
|
||||
real*8 lattice_omega
|
||||
external lattice_omega
|
||||
|
||||
* **** Initializationsr ****
|
||||
call D3dB_n2ft3d(1,n2ft3d)
|
||||
call D3dB_nx(1,nx)
|
||||
call D3dB_ny(1,ny)
|
||||
call D3dB_nz(1,nz)
|
||||
dv = lattice_omega()/dble(nx*ny*nz)
|
||||
|
||||
|
||||
* **** Push memory ****
|
||||
value = MA_push_get(mt_dbl,(3*n2ft3d),'r_grid',
|
||||
> r_grid(2),r_grid(1))
|
||||
value = value.and.
|
||||
> MA_push_get(mt_dbl,(3*n2ft3d),'rho',
|
||||
> rho(2),rho(1))
|
||||
if (.not. value) call errquit(
|
||||
> 'rho_1Fcharge_qmmm: out of stack memory',0, MA_ERR)
|
||||
|
||||
|
||||
* **** Get r_grid and rho ****
|
||||
call lattice_r_grid(dbl_mb(r_grid(1)))
|
||||
call D3dB_rr_Sum(1,dbl_mb(rho1(1)),
|
||||
> dbl_mb(rho1(1)+(ispin-1)*n2ft3d),
|
||||
> dbl_mb(rho(1)))
|
||||
|
||||
* **** Now calculate Fcharge ****
|
||||
call pspw_charge_rho_Fcharge_qmmm(r_charge,q_charge,n_charge,
|
||||
> n2ft3d,
|
||||
> dbl_mb(r_grid(1)),
|
||||
> dbl_mb(rho(1)),dv,
|
||||
> f_charge)
|
||||
|
||||
|
||||
* **** Pop memory ****
|
||||
value = MA_pop_stack(rho(2))
|
||||
value = value.and.MA_pop_stack(r_grid(2))
|
||||
if (.not. value) call errquit(
|
||||
> 'rho_1Fcharge_qmmm: error popping stack memory',0, MA_ERR)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
* ***********************************
|
||||
* * *
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: MOTION.F,v 1.7 2005-10-06 01:04:35 bylaska Exp $
|
||||
* $Id: MOTION.F,v 1.8 2007-03-22 20:46:22 bylaska Exp $
|
||||
|
||||
* Note - The format of the MOTION file was changed on 5-8-2002.
|
||||
*
|
||||
|
|
@ -44,9 +44,9 @@
|
|||
|
||||
* **** external functions ***
|
||||
real*8 lattice_omega
|
||||
integer ion_nion,control_it_out,pspw_qmmm_nion
|
||||
integer ion_nion,control_it_out
|
||||
external lattice_omega
|
||||
external ion_nion,control_it_out,pspw_qmmm_nion
|
||||
external ion_nion,control_it_out
|
||||
|
||||
|
||||
call Parallel_taskid(taskid)
|
||||
|
|
@ -79,10 +79,6 @@
|
|||
> status='new')
|
||||
end if
|
||||
|
||||
c ! taken out on 5-8-2002
|
||||
c open(unit=19,file=full_filename,form='formatted')
|
||||
c write(19,*) control_it_out(),(ion_nion()+pspw_qmmm_nion()),
|
||||
c > lattice_omega()
|
||||
end if
|
||||
|
||||
return
|
||||
|
|
@ -126,20 +122,16 @@ c > lattice_omega()
|
|||
integer taskid,i
|
||||
|
||||
* **** external functions ***
|
||||
logical pspw_qmmm_found
|
||||
integer ion_nion,pspw_qmmm_nion
|
||||
integer ion_nion
|
||||
real*8 ion_rion, ion_vion,lattice_omega
|
||||
real*8 pspw_qmmm_rion,pspw_qmmm_vion
|
||||
external pspw_qmmm_found
|
||||
external ion_nion,pspw_qmmm_nion
|
||||
external ion_nion
|
||||
external ion_rion, ion_vion,lattice_omega
|
||||
external pspw_qmmm_rion,pspw_qmmm_vion
|
||||
|
||||
call Parallel_taskid(taskid)
|
||||
|
||||
if (taskid.eq.MASTER) then
|
||||
write(19,110) time,
|
||||
> (ion_nion()+pspw_qmmm_nion()),
|
||||
> ion_nion(),
|
||||
> lattice_omega()
|
||||
do i=1,ion_nion()
|
||||
write(19,111) ion_rion(1,i),
|
||||
|
|
@ -149,17 +141,6 @@ c > lattice_omega()
|
|||
> ion_vion(2,i),
|
||||
> ion_vion(3,i)
|
||||
end do
|
||||
if (pspw_qmmm_found()) then
|
||||
do i=1,pspw_qmmm_nion()
|
||||
write(19,111) pspw_qmmm_rion(1,i),
|
||||
> pspw_qmmm_rion(2,i),
|
||||
> pspw_qmmm_rion(3,i),
|
||||
> pspw_qmmm_vion(1,i),
|
||||
> pspw_qmmm_vion(2,i),
|
||||
> pspw_qmmm_vion(3,i)
|
||||
|
||||
end do
|
||||
end if
|
||||
call util_flush(19)
|
||||
end if
|
||||
110 format(e14.6,I6,e14.6)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
*
|
||||
* $Id: psp.F,v 1.24 2007-03-19 19:05:38 bylaska Exp $
|
||||
* $Id: psp.F,v 1.25 2007-03-22 20:46:22 bylaska Exp $
|
||||
*
|
||||
|
||||
#define TCGMSG
|
||||
|
|
@ -196,6 +196,21 @@ c value = value.and.MA_free_heap(vnlnrm(2))
|
|||
end
|
||||
|
||||
|
||||
* ***********************************
|
||||
* * *
|
||||
* * psp_zv_ptr *
|
||||
* * *
|
||||
* ***********************************
|
||||
integer function psp_zv_ptr()
|
||||
implicit none
|
||||
|
||||
#include "psp.fh"
|
||||
|
||||
psp_zv_ptr = zv(1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***********************************
|
||||
* * *
|
||||
* * psp_amass *
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
c
|
||||
c $Id: pspw_LJ.F,v 1.10 2003-10-17 22:57:49 carlfahl Exp $
|
||||
c $Id: pspw_LJ.F,v 1.11 2007-03-22 20:46:23 bylaska Exp $
|
||||
c
|
||||
|
||||
* *************************
|
||||
|
|
@ -9,15 +9,12 @@ c
|
|||
* *************************
|
||||
subroutine pspw_LJ_init(rtdb)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
integer rtdb
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "rtdb.fh"
|
||||
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_LJ.fh"
|
||||
#include "errquit.fh"
|
||||
|
||||
* **** local variables
|
||||
logical value
|
||||
|
|
@ -26,61 +23,42 @@ c
|
|||
parameter(MASTER=0)
|
||||
|
||||
integer i,j
|
||||
integer nkatm_mm,nkatm_ion
|
||||
integer nkatm
|
||||
real*8 sigma_epsilon(2)
|
||||
character*30 rtdbname
|
||||
|
||||
* **** external functions ****
|
||||
character*2 pspw_qmmm_atom,ion_atom
|
||||
integer pspw_qmmm_nkatm,ion_nkatm
|
||||
integer pspw_qmmm_katm,ion_katm
|
||||
external pspw_qmmm_atom,ion_atom
|
||||
external pspw_qmmm_nkatm,ion_nkatm
|
||||
external pspw_qmmm_katm,ion_katm
|
||||
character*2 ion_atom
|
||||
integer ion_nkatm
|
||||
integer ion_katm
|
||||
external ion_atom
|
||||
external ion_nkatm
|
||||
external ion_katm
|
||||
|
||||
|
||||
nkatm_mm = pspw_qmmm_nkatm()
|
||||
nkatm_ion = ion_nkatm()
|
||||
nkatm = ion_nkatm()
|
||||
|
||||
* **** allocate LJ parameters ****
|
||||
value = rtdb_parallel(.true.)
|
||||
value = MA_alloc_get(mt_dbl,nkatm,
|
||||
> 'sigma_all',sigma_all(2),sigma_all(1))
|
||||
value = value.and.
|
||||
> MA_alloc_get(mt_dbl,nkatm_mm,
|
||||
> 'sigma_mm',sigma_mm(2),sigma_mm(1))
|
||||
value = value.and.
|
||||
> MA_alloc_get(mt_dbl,nkatm_mm,
|
||||
> 'epsilon_mm',epsilon_mm(2),epsilon_mm(1))
|
||||
value = MA_alloc_get(mt_dbl,nkatm_ion,
|
||||
> 'sigma_ion',sigma_ion(2),sigma_ion(1))
|
||||
value = value.and.
|
||||
> MA_alloc_get(mt_dbl,nkatm_ion,
|
||||
> 'epsilon_ion',epsilon_ion(2),epsilon_ion(1))
|
||||
> MA_alloc_get(mt_dbl,nkatm,
|
||||
> 'epsilon_ion',epsilon_all(2),epsilon_all(1))
|
||||
if (.not. value) call errquit('out of heap memory',0, MA_ERR)
|
||||
|
||||
|
||||
* **** read in MM LJ parameters ***
|
||||
do i=1,nkatm_mm
|
||||
rtdbname = 'pspw_LJ_param_mm:'//pspw_qmmm_atom(i)
|
||||
j = index(rtdbname,' ')-1
|
||||
value = rtdb_get(rtdb,rtdbname(1:j),
|
||||
> mt_dbl,2,sigma_epsilon)
|
||||
if (.not. value) call dcopy(2,0.0d0,0,sigma_epsilon,1)
|
||||
|
||||
dbl_mb(sigma_mm(1)+i-1) =sigma_epsilon(1)/(0.529177d0)
|
||||
dbl_mb(epsilon_mm(1)+i-1)=sigma_epsilon(2)/(23.06d0*27.2116d0)
|
||||
end do
|
||||
|
||||
|
||||
* **** read in Ion LJ parameters ***
|
||||
do i=1,nkatm_ion
|
||||
* **** read in LJ parameters ***
|
||||
do i=1,nkatm
|
||||
rtdbname = 'pspw_LJ_param_ion:'//ion_atom(i)
|
||||
j = index(rtdbname,' ')-1
|
||||
value = rtdb_get(rtdb,rtdbname(1:j),
|
||||
> mt_dbl,2,sigma_epsilon)
|
||||
if (.not. value) call dcopy(2,0.0d0,0,sigma_epsilon,1)
|
||||
|
||||
dbl_mb(sigma_ion(1)+i-1) =sigma_epsilon(1)/(0.529177d0)
|
||||
dbl_mb(epsilon_ion(1)+i-1)=sigma_epsilon(2)/(23.06d0*27.2116d0)
|
||||
dbl_mb(sigma_all(1)+i-1) =sigma_epsilon(1)/(0.529177d0)
|
||||
dbl_mb(epsilon_all(1)+i-1)=sigma_epsilon(2)/(23.06d0*27.2116d0)
|
||||
end do
|
||||
|
||||
|
||||
|
|
@ -88,30 +66,17 @@ c
|
|||
if (taskid.eq.MASTER) then
|
||||
write(*,*)
|
||||
write(*,*) 'Ion LJ Parameters:'
|
||||
do i=1,nkatm_ion
|
||||
do i=1,nkatm
|
||||
write(*,'(A3,4x,A12,E10.3,A9,A12,E10.3,A9)')
|
||||
> ion_atom(i),
|
||||
> 'sigma:',
|
||||
> dbl_mb(sigma_ion(1)+i-1)*(0.529177d0),
|
||||
> dbl_mb(sigma_all(1)+i-1)*(0.529177d0),
|
||||
> 'Angstrom',
|
||||
> 'epsilon:',
|
||||
> dbl_mb(epsilon_ion(1)+ i-1)*(23.06d0*27.2116d0),
|
||||
> dbl_mb(epsilon_all(1)+ i-1)*(23.06d0*27.2116d0),
|
||||
> 'kcal/mol'
|
||||
end do
|
||||
write(*,*)
|
||||
write(*,*) 'MM LJ Parameters:'
|
||||
do i=1,nkatm_mm
|
||||
write(*,'(A3,4x,A12,E10.3,A9,A12,E10.3,A9)')
|
||||
> pspw_qmmm_atom(i),
|
||||
> 'sigma:',
|
||||
> dbl_mb(sigma_mm(1)+i-1)*(0.529177d0),
|
||||
> 'Angstrom',
|
||||
> 'epsilon:',
|
||||
> dbl_mb(epsilon_mm(1)+i-1)*(23.06d0*27.2116d0),
|
||||
> 'kcal/mol'
|
||||
end do
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
return
|
||||
|
|
@ -119,34 +84,30 @@ c
|
|||
|
||||
* *************************
|
||||
* * *
|
||||
* * pspw_LJ_end *
|
||||
* * pspw_LJ_end *
|
||||
* * *
|
||||
* *************************
|
||||
subroutine pspw_LJ_end()
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_LJ.fh"
|
||||
#include "errquit.fh"
|
||||
|
||||
logical value
|
||||
|
||||
value = MA_free_heap(sigma_mm(2))
|
||||
value = value.and.MA_free_heap(epsilon_mm(2))
|
||||
value = value.and.MA_free_heap(sigma_ion(2))
|
||||
value = value.and.MA_free_heap(epsilon_ion(2))
|
||||
if (.not.value) call errquit('error MA_free_heap',0, MA_ERR)
|
||||
|
||||
value = MA_free_heap(sigma_all(2))
|
||||
value = value.and.MA_free_heap(epsilon_all(2))
|
||||
if (.not.value) call errquit('pspw_LJ_end: error MA_free_heap',
|
||||
> 0,MA_ERR)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * *
|
||||
* * pspw_LJ_param_Ion *
|
||||
* * *
|
||||
* * *
|
||||
* ***************************
|
||||
subroutine pspw_LJ_param_Ion(i,s,e)
|
||||
implicit none
|
||||
|
|
@ -154,8 +115,6 @@ c
|
|||
real*8 s,e
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_LJ.fh"
|
||||
|
||||
* **** local variables ****
|
||||
|
|
@ -166,357 +125,253 @@ c
|
|||
external ion_katm
|
||||
|
||||
ka = ion_katm(i)
|
||||
s = dbl_mb(sigma_ion(1) + ka-1)
|
||||
e = dbl_mb(epsilon_ion(1)+ ka-1)
|
||||
s = dbl_mb(sigma_all(1) + ka-1)
|
||||
e = dbl_mb(epsilon_all(1)+ ka-1)
|
||||
return
|
||||
end
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_LJ_param_MM *
|
||||
* * *
|
||||
* ***************************
|
||||
subroutine pspw_LJ_param_MM(i,s,e)
|
||||
|
||||
c *************************************
|
||||
c * *
|
||||
c * pspw_LJ_E *
|
||||
c * *
|
||||
c *************************************
|
||||
real*8 function pspw_LJ_E(nion,nion_qm,katm,
|
||||
> nfrag,indx_frag_start,size_frag,
|
||||
> nshl3d,rcell,
|
||||
> rion)
|
||||
implicit none
|
||||
integer i
|
||||
real*8 s,e
|
||||
integer nion,nion_qm
|
||||
integer katm(*)
|
||||
integer nfrag
|
||||
integer indx_frag_start(*),size_frag(*)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
real*8 rion(3,*)
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_LJ.fh"
|
||||
|
||||
* **** local variables ****
|
||||
integer ka
|
||||
|
||||
* **** external functions ****
|
||||
integer pspw_qmmm_katm
|
||||
external pspw_qmmm_katm
|
||||
|
||||
ka = pspw_qmmm_katm(i)
|
||||
s = dbl_mb(sigma_mm(1) + ka-1)
|
||||
e = dbl_mb(epsilon_mm(1)+ ka-1)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_LJ_Energy_MM *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function pspw_LJ_Energy_MM()
|
||||
implicit none
|
||||
|
||||
* **** local variables ****
|
||||
real*8 E,si,ei,sj,ej
|
||||
real*8 ri(3),rj(3)
|
||||
integer i,j,nion
|
||||
|
||||
* **** external functions ****
|
||||
logical pspw_qmmm_different_frag
|
||||
integer pspw_qmmm_nion
|
||||
real*8 pspw_qmmm_rion,LJ_Repulsion
|
||||
external pspw_qmmm_different_frag
|
||||
external pspw_qmmm_nion
|
||||
external pspw_qmmm_rion,LJ_Repulsion
|
||||
|
||||
nion = pspw_qmmm_nion()
|
||||
|
||||
E = 0.0d0
|
||||
do i=1,nion
|
||||
call pspw_LJ_param_MM(i,si,ei)
|
||||
if (ei.gt.1.0d-10) then
|
||||
ri(1) = pspw_qmmm_rion(1,i)
|
||||
ri(2) = pspw_qmmm_rion(2,i)
|
||||
ri(3) = pspw_qmmm_rion(3,i)
|
||||
do j=i+1,nion
|
||||
if (pspw_qmmm_different_frag(i,j)) then
|
||||
call pspw_LJ_param_MM(j,sj,ej)
|
||||
if (ej.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
E = E + LJ_Repulsion(ri,si,ei,rj,sj,ej)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
|
||||
pspw_LJ_Energy_MM = E
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_LJ_Energy_IonMM *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function pspw_LJ_Energy_IonMM()
|
||||
implicit none
|
||||
|
||||
* **** local variables ****
|
||||
real*8 E,si,ei,sj,ej
|
||||
real*8 ri(3),rj(3)
|
||||
integer i,j,nion_ion,nion_mm
|
||||
|
||||
* **** external functions ****
|
||||
integer pspw_qmmm_nion,ion_nion
|
||||
real*8 pspw_qmmm_rion,ion_rion,LJ_Repulsion
|
||||
external pspw_qmmm_nion,ion_nion
|
||||
external pspw_qmmm_rion,ion_rion,LJ_Repulsion
|
||||
|
||||
nion_mm = pspw_qmmm_nion()
|
||||
nion_ion = ion_nion()
|
||||
|
||||
E = 0.0d0
|
||||
do i=1,nion_ion
|
||||
call pspw_LJ_param_Ion(i,si,ei)
|
||||
if (ei.gt.1.0d-10) then
|
||||
ri(1) = ion_rion(1,i)
|
||||
ri(2) = ion_rion(2,i)
|
||||
ri(3) = ion_rion(3,i)
|
||||
do j=1,nion_mm
|
||||
call pspw_LJ_param_MM(j,sj,ej)
|
||||
if (ej.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
E = E + LJ_Repulsion(ri,si,ei,rj,sj,ej)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
|
||||
pspw_LJ_Energy_IonMM = E
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_LJ_Fmm *
|
||||
* * *
|
||||
* ***************************
|
||||
|
||||
* This routine calculates the force on the MM ions from the
|
||||
* Lenard-Jones interaction with the MM atoms and the QM ions and
|
||||
* then adds it to the MM ion force, f_mm.
|
||||
*
|
||||
* Exit - f_mm: force on the MM ions
|
||||
*
|
||||
* Uses - ion_nion,ion_rion
|
||||
* pspw_qmmm_nion,pspw_qmmm_rion
|
||||
* pspw_qmmm_different_frag
|
||||
*
|
||||
* Author - Eric Bylaska
|
||||
*
|
||||
subroutine pspw_LJ_Fmm(f_mm)
|
||||
implicit none
|
||||
real*8 f_mm(3,*)
|
||||
|
||||
* **** local variables ****
|
||||
real*8 si,ei,sj,ej
|
||||
real*8 ri(3),rj(3),fion(3)
|
||||
integer i,j
|
||||
integer w1,a,k1,kk1
|
||||
integer w2,b,k2,kk2
|
||||
real*8 e1,s1,e2,s2,E
|
||||
|
||||
* **** external functions ****
|
||||
logical pspw_qmmm_different_frag
|
||||
integer pspw_qmmm_nion,ion_nion
|
||||
real*8 pspw_qmmm_rion,ion_rion
|
||||
external pspw_qmmm_different_frag
|
||||
external pspw_qmmm_nion,ion_nion
|
||||
external pspw_qmmm_rion,ion_rion
|
||||
real*8 pspw_VLJ_E_periodic
|
||||
external pspw_VLJ_E_periodic
|
||||
|
||||
|
||||
E = 0.0d0
|
||||
|
||||
* **** force from the MM atoms ****
|
||||
do i=1,pspw_qmmm_nion()
|
||||
call pspw_LJ_param_MM(i,si,ei)
|
||||
if (ei.gt.1.0d-10) then
|
||||
ri(1) = pspw_qmmm_rion(1,i)
|
||||
ri(2) = pspw_qmmm_rion(2,i)
|
||||
ri(3) = pspw_qmmm_rion(3,i)
|
||||
do j=i+1,pspw_qmmm_nion()
|
||||
if (pspw_qmmm_different_frag(i,j)) then
|
||||
call pspw_LJ_param_MM(j,sj,ej)
|
||||
if (ej.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
call LJ_Ion_Repulsion_Force(ri,si,ei,f_mm(1,i),
|
||||
> rj,sj,ej,f_mm(1,j))
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
|
||||
* **** force from the QM atoms ****
|
||||
call dcopy(3,0.0d0,0,fion,1)
|
||||
do i=1,ion_nion()
|
||||
call pspw_LJ_param_Ion(i,si,ei)
|
||||
if (ei.gt.1.0d-10) then
|
||||
ri(1) = ion_rion(1,i)
|
||||
ri(2) = ion_rion(2,i)
|
||||
ri(3) = ion_rion(3,i)
|
||||
do j=1,pspw_qmmm_nion()
|
||||
call pspw_LJ_param_MM(j,sj,ej)
|
||||
if (ej.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
call LJ_Ion_Repulsion_Force(ri,si,ei,fion,
|
||||
> rj,sj,ej,f_mm(1,j))
|
||||
c **** QM/MM LJ energy ****
|
||||
do i=1,nion_qm
|
||||
e1 = dbl_mb(epsilon_all(1)+katm(i)-1)
|
||||
s1 = dbl_mb(sigma_all(1) +katm(i)-1)
|
||||
do j = nion_qm+1,nion
|
||||
e2 = dbl_mb(epsilon_all(1)+katm(j)-1)
|
||||
if ((e1.gt.1.0d-11).and.(e2.gt.1.0d-11)) then
|
||||
s2 = dbl_mb(sigma_all(1)+katm(j)-1)
|
||||
E = E + pspw_VLJ_E_periodic(e1,s1,rion(1,i),
|
||||
> e2,s2,rion(1,j),
|
||||
> nshl3d,rcell)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
c **** MM/MM LJ energy ****
|
||||
do w1 = 1,nfrag-1
|
||||
do w2 = w1+1,nfrag
|
||||
k1 = indx_frag_start(w1)
|
||||
k2 = indx_frag_start(w2)
|
||||
kk1 = k1
|
||||
do a=1,size_frag(w1)
|
||||
kk2 = k2
|
||||
do b=1,size_frag(w2)
|
||||
e1 = dbl_mb(epsilon_all(1)+katm(kk1)-1)
|
||||
e2 = dbl_mb(epsilon_all(1)+katm(kk2)-1)
|
||||
if ((e1.gt.1.0d-11).and.(e2.gt.1.0d-11)) then
|
||||
s1 = dbl_mb(sigma_all(1)+katm(kk1)-1)
|
||||
s2 = dbl_mb(sigma_all(1)+katm(kk2)-1)
|
||||
E = E + pspw_VLJ_E_periodic(e1,s1,rion(1,kk1),
|
||||
> e2,s2,rion(1,kk2),
|
||||
> nshl3d,rcell)
|
||||
end if
|
||||
kk2 = kk2 + 1
|
||||
end do
|
||||
kk1 = kk1 + 1
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
pspw_LJ_E = E
|
||||
return
|
||||
end
|
||||
|
||||
c *********************************************
|
||||
c * *
|
||||
c * pspw_VLJ_E_periodic *
|
||||
c * *
|
||||
c *********************************************
|
||||
|
||||
real*8 function pspw_VLJ_E_periodic(e1,s1,r1,e2,s2,r2,
|
||||
> nshl3d,rcell)
|
||||
implicit none
|
||||
real*8 e1,s1,r1(3)
|
||||
real*8 e2,s2,r2(3)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
|
||||
* **** local variables ****
|
||||
integer l
|
||||
real*8 dx,dy,dz
|
||||
real*8 x,y,z,r
|
||||
real*8 sigma12,epsilon12
|
||||
real*8 u,u6,u12,E
|
||||
|
||||
E = 0.0d0
|
||||
sigma12 = dsqrt(s1*s2)
|
||||
epsilon12 = dsqrt(e1*e2)
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
do l=1,nshl3d
|
||||
x = dx + rcell(l,1)
|
||||
y = dy + rcell(l,2)
|
||||
z = dz + rcell(l,3)
|
||||
r = dsqrt(x**2 + y**2 + z**2)
|
||||
u = (sigma12/r)
|
||||
u6 = u**6
|
||||
u12 = u6**2
|
||||
E = E + 4.0d0*epsilon12*(u12-u6)
|
||||
end do
|
||||
|
||||
pspw_VLJ_E_periodic = E
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
c *************************************
|
||||
c * *
|
||||
c * pspw_LJ_fion *
|
||||
c * *
|
||||
c *************************************
|
||||
subroutine pspw_LJ_fion(nion,nion_qm,katm,
|
||||
> nfrag,indx_frag_start,size_frag,
|
||||
> nshl3d,rcell,
|
||||
> rion,fion)
|
||||
implicit none
|
||||
integer nion,nion_qm
|
||||
integer katm(*)
|
||||
integer nfrag
|
||||
integer indx_frag_start(*),size_frag(*)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
real*8 rion(3,*)
|
||||
real*8 fion(3,*)
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "pspw_LJ.fh"
|
||||
|
||||
* **** local variables ****
|
||||
integer i,j
|
||||
integer w1,a,k1,kk1
|
||||
integer w2,b,k2,kk2
|
||||
real*8 e1,s1,e2,s2
|
||||
|
||||
|
||||
c **** QM/MM LJ energy ****
|
||||
do i=1,nion_qm
|
||||
e1 = dbl_mb(epsilon_all(1)+katm(i)-1)
|
||||
s1 = dbl_mb(sigma_all(1) +katm(i)-1)
|
||||
do j = nion_qm+1,nion
|
||||
e2 = dbl_mb(epsilon_all(1)+katm(j)-1)
|
||||
if ((e1.gt.1.0d-11).and.(e2.gt.1.0d-11)) then
|
||||
s2 = dbl_mb(sigma_all(1)+katm(j)-1)
|
||||
call pspw_VLJ_fion_periodic(e1,s1,rion(1,i),fion(1,i),
|
||||
> e2,s2,rion(1,j),fion(1,j),
|
||||
> nshl3d,rcell)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
c **** MM/MM LJ energy ****
|
||||
do w1 = 1,nfrag-1
|
||||
do w2 = w1+1,nfrag
|
||||
k1 = indx_frag_start(w1)
|
||||
k2 = indx_frag_start(w2)
|
||||
kk1 = k1
|
||||
do a=1,size_frag(w1)
|
||||
kk2 = k2
|
||||
do b=1,size_frag(w2)
|
||||
e1 = dbl_mb(epsilon_all(1)+katm(kk1)-1)
|
||||
e2 = dbl_mb(epsilon_all(1)+katm(kk2)-1)
|
||||
if ((e1.gt.1.0d-11).and.(e2.gt.1.0d-11)) then
|
||||
s1 = dbl_mb(sigma_all(1)+katm(kk1)-1)
|
||||
s2 = dbl_mb(sigma_all(1)+katm(kk2)-1)
|
||||
call pspw_VLJ_fion_periodic(e1,s1,rion(1,kk1),fion(1,kk1),
|
||||
> e2,s2,rion(1,kk2),fion(1,kk2),
|
||||
> nshl3d,rcell)
|
||||
end if
|
||||
kk2 = kk2 + 1
|
||||
end do
|
||||
kk1 = kk1 + 1
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_LJ_Fion *
|
||||
* * *
|
||||
* ***************************
|
||||
c *********************************************
|
||||
c * *
|
||||
c * pspw_VLJ_fion_periodic *
|
||||
c * *
|
||||
c *********************************************
|
||||
|
||||
* This routine calculates the force on the QM ions from the
|
||||
* Lenard-Jones interaction with the MM atoms and adds it
|
||||
* to the QM ion force, fion.
|
||||
*
|
||||
* Exit - fion: force on the QM ions
|
||||
*
|
||||
* Uses - ion_nion,ion_rion
|
||||
* pspw_qmmm_nion,pspw_qmmm_rion
|
||||
*
|
||||
* Author - Eric Bylaska
|
||||
|
||||
subroutine pspw_LJ_Fion(fion)
|
||||
subroutine pspw_VLJ_fion_periodic(e1,s1,r1,f1,
|
||||
> e2,s2,r2,f2,
|
||||
> nshl3d,rcell)
|
||||
implicit none
|
||||
real*8 fion(3,*)
|
||||
real*8 e1,s1,r1(3),f1(3)
|
||||
real*8 e2,s2,r2(3),f2(3)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
|
||||
* **** local variables ****
|
||||
real*8 si,ei,sj,ej
|
||||
real*8 ri(3),rj(3),f_mm(3)
|
||||
integer i,j,nion_ion,nion_mm
|
||||
integer l
|
||||
real*8 dx,dy,dz
|
||||
real*8 x,y,z,r
|
||||
real*8 sigma12,epsilon12
|
||||
real*8 u,u6,u12,dVLJ
|
||||
|
||||
* **** external functions ****
|
||||
integer pspw_qmmm_nion,ion_nion
|
||||
real*8 pspw_qmmm_rion,ion_rion
|
||||
external pspw_qmmm_nion,ion_nion
|
||||
external pspw_qmmm_rion,ion_rion
|
||||
|
||||
nion_mm = pspw_qmmm_nion()
|
||||
nion_ion = ion_nion()
|
||||
|
||||
call dcopy(3,0.0d0,0,f_mm,1)
|
||||
do i=1,nion_ion
|
||||
call pspw_LJ_param_Ion(i,si,ei)
|
||||
if (ei.gt.1.0d-10) then
|
||||
ri(1) = ion_rion(1,i)
|
||||
ri(2) = ion_rion(2,i)
|
||||
ri(3) = ion_rion(3,i)
|
||||
do j=1,nion_mm
|
||||
call pspw_LJ_param_MM(j,sj,ej)
|
||||
if (ej.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
call LJ_Ion_Repulsion_Force(ri,si,ei,fion(1,i),
|
||||
> rj,sj,ej,f_mm)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
sigma12 = dsqrt(s1*s2)
|
||||
epsilon12 = dsqrt(e1*e2)
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
do l=1,nshl3d
|
||||
x = dx + rcell(l,1)
|
||||
y = dy + rcell(l,2)
|
||||
z = dz + rcell(l,3)
|
||||
r = dsqrt(x**2 + y**2 + z**2)
|
||||
u = (sigma12/r)
|
||||
u6 = u**6
|
||||
u12 = u6**2
|
||||
dVLJ = -(4.0d0*epsilon12/r)*(12.0d0*u12-6.0d0*u6)
|
||||
f1(1) = f1(1) - (x/r)*dVLJ
|
||||
f1(2) = f1(2) - (y/r)*dVLJ
|
||||
f1(3) = f1(3) - (z/r)*dVLJ
|
||||
f2(1) = f2(1) + (x/r)*dVLJ
|
||||
f2(2) = f2(2) + (y/r)*dVLJ
|
||||
f2(3) = f2(3) + (z/r)*dVLJ
|
||||
end do
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* **********************************
|
||||
* * *
|
||||
* * LJ_Repulsion *
|
||||
* * *
|
||||
* **********************************
|
||||
|
||||
* This routine calculates the LJ Repulsion between two centers
|
||||
*
|
||||
real*8 function LJ_Repulsion(r1,s1,e1,r2,s2,e2)
|
||||
implicit none
|
||||
real*8 r1(3),s1,e1
|
||||
real*8 r2(3),s2,e2
|
||||
|
||||
* **** local variables ****
|
||||
real*8 xo,yo,zo,d2
|
||||
real*8 ee,ss
|
||||
|
||||
ee = dsqrt(e1*e2)
|
||||
ss = (s1*s2)
|
||||
|
||||
xo = r2(1)-r1(1)
|
||||
yo = r2(2)-r1(2)
|
||||
zo = r2(3)-r1(3)
|
||||
|
||||
d2 = (xo*xo + yo*yo + zo*zo)
|
||||
|
||||
|
||||
LJ_Repulsion = 4.0d0*ee * ((ss/d2)**6 - (ss/d2)**3)
|
||||
return
|
||||
end
|
||||
|
||||
* **********************************
|
||||
* * *
|
||||
* * LJ_Repulsion_Force *
|
||||
* * *
|
||||
* **********************************
|
||||
|
||||
* This routine calculates the LJ Repulsion force between two centers
|
||||
*
|
||||
subroutine LJ_Ion_Repulsion_Force(r1,s1,e1,f1,r2,s2,e2,f2)
|
||||
implicit none
|
||||
real*8 r1(3),s1,e1
|
||||
real*8 f1(3)
|
||||
real*8 r2(3),s2,e2
|
||||
real*8 f2(3)
|
||||
|
||||
* **** local variables ****
|
||||
real*8 xo,yo,zo,d2
|
||||
real*8 der
|
||||
real*8 ee,ss,fx,fy,fz
|
||||
|
||||
ee = dsqrt(e1*e2)
|
||||
ss = (s1*s2)
|
||||
|
||||
xo = r2(1)-r1(1)
|
||||
yo = r2(2)-r1(2)
|
||||
zo = r2(3)-r1(3)
|
||||
|
||||
d2 = (xo*xo + yo*yo + zo*zo)
|
||||
|
||||
der = (4.0d0*ee/d2) * ( -12.0d0*(ss/d2)**6
|
||||
> + 6.0d0*(ss/d2)**3 )
|
||||
fx = -(xo)*der
|
||||
fy = -(yo)*der
|
||||
fz = -(zo)*der
|
||||
|
||||
f2(1) = f2(1) + fx
|
||||
f2(2) = f2(2) + fy
|
||||
f2(3) = f2(3) + fz
|
||||
|
||||
f1(1) = f1(1) - fx
|
||||
f1(2) = f1(2) - fy
|
||||
f1(3) = f1(3) - fz
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,12 +1,9 @@
|
|||
c
|
||||
c $Id: pspw_LJ.fh,v 1.4 2001-11-01 16:37:12 bylaska Exp $
|
||||
c $Id: pspw_LJ.fh,v 1.5 2007-03-22 20:46:23 bylaska Exp $
|
||||
c
|
||||
|
||||
* ***** pspw_LJ common block ************************************
|
||||
integer sigma_mm(2),epsilon_mm(2)
|
||||
integer sigma_ion(2),epsilon_ion(2)
|
||||
common / pspw_LJ / sigma_mm, ! MM LJ parameters
|
||||
> epsilon_mm, ! MM LJ parameters
|
||||
> sigma_ion, ! Ion LJ parameters
|
||||
> epsilon_ion ! Ion LJ parameters
|
||||
integer sigma_all(2),epsilon_all(2)
|
||||
common / pspw_LJ / sigma_all, ! LJ parameters
|
||||
> epsilon_all ! LJ parameters
|
||||
* ******************************************************************
|
||||
|
|
|
|||
|
|
@ -1,441 +1,9 @@
|
|||
c
|
||||
c $Id: pspw_Q.F,v 1.11 2003-10-17 22:57:53 carlfahl Exp $
|
||||
c $Id: pspw_Q.F,v 1.12 2007-03-22 20:46:23 bylaska Exp $
|
||||
c
|
||||
|
||||
|
||||
* *************************
|
||||
* * *
|
||||
* * pspw_Q_init *
|
||||
* * *
|
||||
* *************************
|
||||
subroutine pspw_Q_init(rtdb)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
integer rtdb
|
||||
|
||||
#include "mafdecls.fh"
|
||||
#include "rtdb.fh"
|
||||
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_Q.fh"
|
||||
|
||||
* **** local variables
|
||||
logical value
|
||||
integer taskid
|
||||
integer MASTER
|
||||
parameter(MASTER=0)
|
||||
|
||||
integer i,j
|
||||
integer nkatm_mm,nkatm_ion
|
||||
real*8 q
|
||||
character*30 rtdbname
|
||||
|
||||
* **** external functions ****
|
||||
character*2 pspw_qmmm_atom,ion_atom
|
||||
integer pspw_qmmm_nkatm,ion_nkatm
|
||||
integer pspw_qmmm_katm,ion_katm
|
||||
real*8 psp_zv
|
||||
external pspw_qmmm_atom,ion_atom
|
||||
external pspw_qmmm_nkatm,ion_nkatm
|
||||
external pspw_qmmm_katm,ion_katm
|
||||
external psp_zv
|
||||
|
||||
|
||||
nkatm_mm = pspw_qmmm_nkatm()
|
||||
nkatm_ion = ion_nkatm()
|
||||
|
||||
|
||||
* **** allocate Electrostatic parameters ****
|
||||
* **** default is zero ****
|
||||
value = rtdb_parallel(.true.)
|
||||
value = value.and.
|
||||
> MA_alloc_get(mt_dbl,nkatm_mm,
|
||||
> 'q_mm',q_mm(2),q_mm(1))
|
||||
value = MA_alloc_get(mt_dbl,nkatm_ion,
|
||||
> 'q_ion',q_ion(2),q_ion(1))
|
||||
if (.not. value) call errquit('out of heap memory',0, MA_ERR)
|
||||
|
||||
|
||||
* **** read in MM Electrostatic parameters ***
|
||||
do i=1,nkatm_mm
|
||||
rtdbname = 'pspw_Q_param_mm:'//pspw_qmmm_atom(i)
|
||||
j = index(rtdbname,' ')-1
|
||||
value = rtdb_get(rtdb,rtdbname(1:j),
|
||||
> mt_dbl,1,q)
|
||||
if (.not. value) q=0.0d0
|
||||
|
||||
dbl_mb(q_mm(1)+i-1) = q
|
||||
end do
|
||||
|
||||
* **** read in Ion Electrostatic parameters ****
|
||||
* **** default is core charge of pseudopotential ****
|
||||
do i=1,nkatm_ion
|
||||
rtdbname = 'pspw_Q_param_ion:'//ion_atom(i)
|
||||
j = index(rtdbname,' ')-1
|
||||
value = rtdb_get(rtdb,rtdbname(1:j),
|
||||
> mt_dbl,1,q)
|
||||
if (.not. value) q=psp_zv(i)
|
||||
|
||||
dbl_mb(q_ion(1)+i-1) = q
|
||||
end do
|
||||
|
||||
call Parallel_taskid(taskid)
|
||||
if (taskid.eq.MASTER) then
|
||||
write(*,*)
|
||||
write(*,*) 'Ion Electrostatic Parameters:'
|
||||
do i=1,nkatm_ion
|
||||
write(*,'(A3,4x,A12,E10.3,A3)')
|
||||
> ion_atom(i),
|
||||
> 'Q:',dbl_mb(q_ion(1) + i-1),'au'
|
||||
end do
|
||||
write(*,*)
|
||||
write(*,*) 'MM Electrostatic Parameters:'
|
||||
do i=1,nkatm_mm
|
||||
write(*,'(A3,4x,A12,E10.3,A3)')
|
||||
> pspw_qmmm_atom(i),
|
||||
> 'Q:',dbl_mb(q_mm(1) + i-1),'au'
|
||||
end do
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
* *************************
|
||||
* * *
|
||||
* * pspw_Q_end *
|
||||
* * *
|
||||
* *************************
|
||||
subroutine pspw_Q_end()
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_Q.fh"
|
||||
|
||||
logical value
|
||||
|
||||
value = MA_free_heap(q_mm(2))
|
||||
value = value.and.MA_free_heap(q_ion(2))
|
||||
if (.not.value) call errquit('error MA_free_heap',0, MA_ERR)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_param_Ion *
|
||||
* * *
|
||||
* ***************************
|
||||
subroutine pspw_Q_param_Ion(i,q)
|
||||
implicit none
|
||||
integer i
|
||||
real*8 q
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_Q.fh"
|
||||
|
||||
* **** local variables ****
|
||||
integer ka
|
||||
|
||||
* **** external functions ****
|
||||
integer ion_katm
|
||||
external ion_katm
|
||||
|
||||
ka = ion_katm(i)
|
||||
q = dbl_mb(q_ion(1) + ka-1)
|
||||
return
|
||||
end
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_param_MM *
|
||||
* * *
|
||||
* ***************************
|
||||
subroutine pspw_Q_param_MM(i,q)
|
||||
implicit none
|
||||
integer i
|
||||
real*8 q
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_Q.fh"
|
||||
|
||||
* **** local variables ****
|
||||
integer ka
|
||||
|
||||
* **** external functions ****
|
||||
integer pspw_qmmm_katm
|
||||
external pspw_qmmm_katm
|
||||
|
||||
ka = pspw_qmmm_katm(i)
|
||||
q = dbl_mb(q_mm(1) + ka-1)
|
||||
return
|
||||
end
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_zv_MM *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function pspw_Q_zv_MM(ia)
|
||||
implicit none
|
||||
integer ia
|
||||
|
||||
#include "mafdecls.fh"
|
||||
|
||||
***** LJparam common block ****
|
||||
#include "pspw_Q.fh"
|
||||
|
||||
pspw_Q_zv_MM = dbl_mb(q_mm(1) + ia-1)
|
||||
return
|
||||
end
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_Energy_MM *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function pspw_Q_Energy_MM()
|
||||
implicit none
|
||||
|
||||
* **** local variables ****
|
||||
real*8 E,qi,qj
|
||||
real*8 ri(3),rj(3)
|
||||
integer i,j,nion
|
||||
|
||||
* **** external functions ****
|
||||
logical pspw_qmmm_different_frag
|
||||
integer pspw_qmmm_nion
|
||||
real*8 pspw_qmmm_rion,Q_Electrostatic
|
||||
external pspw_qmmm_different_frag
|
||||
external pspw_qmmm_nion
|
||||
external pspw_qmmm_rion,Q_Electrostatic
|
||||
|
||||
nion = pspw_qmmm_nion()
|
||||
|
||||
E = 0.0d0
|
||||
do i=1,nion
|
||||
call pspw_Q_param_MM(i,qi)
|
||||
if (qi.gt.1.0d-10) then
|
||||
ri(1) = pspw_qmmm_rion(1,i)
|
||||
ri(2) = pspw_qmmm_rion(2,i)
|
||||
ri(3) = pspw_qmmm_rion(3,i)
|
||||
do j=i+1,nion
|
||||
if (pspw_qmmm_different_frag(i,j)) then
|
||||
call pspw_Q_param_MM(j,qj)
|
||||
if (qj.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
E = E + Q_Electrostatic(ri,qi,rj,qj)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
|
||||
pspw_Q_Energy_MM = E
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_Energy_IonMM *
|
||||
* * *
|
||||
* ***************************
|
||||
real*8 function pspw_Q_Energy_IonMM()
|
||||
implicit none
|
||||
|
||||
* **** local variables ****
|
||||
real*8 E,qi,qj
|
||||
real*8 ri(3),rj(3)
|
||||
integer i,j,nion_ion,nion_mm
|
||||
|
||||
* **** external functions ****
|
||||
integer pspw_qmmm_nion,ion_nion
|
||||
real*8 pspw_qmmm_rion,ion_rion,Q_Electrostatic
|
||||
external pspw_qmmm_nion,ion_nion
|
||||
external pspw_qmmm_rion,ion_rion,Q_Electrostatic
|
||||
|
||||
nion_mm = pspw_qmmm_nion()
|
||||
nion_ion = ion_nion()
|
||||
|
||||
E = 0.0d0
|
||||
do i=1,nion_ion
|
||||
call pspw_Q_param_Ion(i,qi)
|
||||
if (qi.gt.1.0d-10) then
|
||||
ri(1) = ion_rion(1,i)
|
||||
ri(2) = ion_rion(2,i)
|
||||
ri(3) = ion_rion(3,i)
|
||||
do j=1,nion_mm
|
||||
call pspw_Q_param_MM(j,qj)
|
||||
if (qj.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
E = E + Q_Electrostatic(ri,qi,rj,qj)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
|
||||
pspw_Q_Energy_IonMM = E
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_Fmm *
|
||||
* * *
|
||||
* ***************************
|
||||
|
||||
* This routine calculates the force on the MM ions from the
|
||||
* Point charge interaction with the MM and QM atoms and adds it
|
||||
* to the MM ion force, f_mm.
|
||||
*
|
||||
* Exit - f_mm: force on the MM ions
|
||||
*
|
||||
* Uses - ion_nion,ion_rion
|
||||
* pspw_qmmm_nion,pspw_qmmm_rion
|
||||
*
|
||||
* Author - Eric Bylaska
|
||||
*
|
||||
subroutine pspw_Q_Fmm(f_mm)
|
||||
implicit none
|
||||
real*8 f_mm(3,*)
|
||||
|
||||
* **** local variables ****
|
||||
real*8 qi,qj
|
||||
real*8 ri(3),rj(3),fion(3)
|
||||
integer i,j
|
||||
|
||||
* **** external functions ****
|
||||
logical pspw_qmmm_different_frag
|
||||
integer pspw_qmmm_nion,ion_nion
|
||||
real*8 pspw_qmmm_rion,ion_rion
|
||||
external pspw_qmmm_different_frag
|
||||
external pspw_qmmm_nion,ion_nion
|
||||
external pspw_qmmm_rion,ion_rion
|
||||
|
||||
|
||||
* **** forces from MM ions ****
|
||||
do i=1,pspw_qmmm_nion()
|
||||
call pspw_Q_param_MM(i,qi)
|
||||
if (qi.gt.1.0d-10) then
|
||||
ri(1) = pspw_qmmm_rion(1,i)
|
||||
ri(2) = pspw_qmmm_rion(2,i)
|
||||
ri(3) = pspw_qmmm_rion(3,i)
|
||||
do j=i+1,pspw_qmmm_nion()
|
||||
if (pspw_qmmm_different_frag(i,j)) then
|
||||
call pspw_Q_param_MM(j,qj)
|
||||
if (qj.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
call Q_Electrostatic_Force(ri,qi,f_mm(1,i),
|
||||
> rj,qj,f_mm(1,j))
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
* **** forces from QM ions ****
|
||||
call dcopy(3,0.0d0,0,fion,1)
|
||||
do i=1,ion_nion()
|
||||
call pspw_Q_param_Ion(i,qi)
|
||||
if (qi.gt.1.0d-10) then
|
||||
ri(1) = ion_rion(1,i)
|
||||
ri(2) = ion_rion(2,i)
|
||||
ri(3) = ion_rion(3,i)
|
||||
do j=1,pspw_qmmm_nion()
|
||||
call pspw_Q_param_MM(j,qj)
|
||||
if (qj.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
call Q_Electrostatic_Force(ri,qi,fion,
|
||||
> rj,qj,f_mm(1,j))
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ***************************
|
||||
* * *
|
||||
* * pspw_Q_Fion *
|
||||
* * *
|
||||
* ***************************
|
||||
|
||||
* This routine calculates the force on the QM ions from the
|
||||
* Point charge interaction with the MM atoms and adds it
|
||||
* to the QM ion force, fion.
|
||||
*
|
||||
* Exit - fion: force on the QM ions
|
||||
*
|
||||
* Uses - ion_nion,ion_rion
|
||||
* pspw_qmmm_nion,pspw_qmmm_rion
|
||||
*
|
||||
* Author - Eric Bylaska
|
||||
*
|
||||
subroutine pspw_Q_Fion(fion)
|
||||
implicit none
|
||||
real*8 fion(3,*)
|
||||
|
||||
* **** local variables ****
|
||||
real*8 qi,qj
|
||||
real*8 ri(3),rj(3),f_mm(3)
|
||||
integer i,j,nion_ion,nion_mm
|
||||
|
||||
* **** external functions ****
|
||||
integer pspw_qmmm_nion,ion_nion
|
||||
real*8 pspw_qmmm_rion,ion_rion
|
||||
external pspw_qmmm_nion,ion_nion
|
||||
external pspw_qmmm_rion,ion_rion
|
||||
|
||||
nion_mm = pspw_qmmm_nion()
|
||||
nion_ion = ion_nion()
|
||||
|
||||
call dcopy(3,0.0d0,0,f_mm,1)
|
||||
do i=1,nion_ion
|
||||
call pspw_Q_param_Ion(i,qi)
|
||||
if (qi.gt.1.0d-10) then
|
||||
ri(1) = ion_rion(1,i)
|
||||
ri(2) = ion_rion(2,i)
|
||||
ri(3) = ion_rion(3,i)
|
||||
do j=1,nion_mm
|
||||
call pspw_Q_param_MM(j,qj)
|
||||
if (qj.gt.1.0d-10) then
|
||||
rj(1) = pspw_qmmm_rion(1,j)
|
||||
rj(2) = pspw_qmmm_rion(2,j)
|
||||
rj(3) = pspw_qmmm_rion(3,j)
|
||||
call Q_Electrostatic_Force(ri,qi,fion(1,i),
|
||||
> rj,qj,f_mm)
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* **********************************
|
||||
|
|
@ -510,3 +78,522 @@ c
|
|||
|
||||
|
||||
|
||||
|
||||
c *************************************
|
||||
c * *
|
||||
c * pspw_Q_E *
|
||||
c * *
|
||||
c *************************************
|
||||
real*8 function pspw_Q_E(nion,nion_qm,katm,psp_zv,amass,
|
||||
> nfrag,indx_frag_start,size_frag,
|
||||
> kfrag,
|
||||
> switch_Rin,switch_Rout,
|
||||
> nshl3d,rcell,
|
||||
> rion)
|
||||
implicit none
|
||||
integer nion,nion_qm
|
||||
integer katm(*)
|
||||
real*8 psp_zv(*),amass(*)
|
||||
integer nfrag
|
||||
integer indx_frag_start(*),size_frag(*)
|
||||
integer kfrag(*)
|
||||
real*8 switch_Rin(*),switch_Rout(*)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
real*8 rion(3,*)
|
||||
|
||||
|
||||
* **** local variables ****
|
||||
integer w1,a,ks1,kk1,n1
|
||||
integer w2,b,ks2,kk2,n2
|
||||
real*8 Rin1,Rout1,rw1_cm(3)
|
||||
real*8 Rin2,Rout2,rw2_cm(3)
|
||||
real*8 Rin,Rout,E
|
||||
|
||||
* **** external functions ****
|
||||
real*8 pspw_Q_E_frag_frag,Q_Electrostatic
|
||||
external pspw_Q_E_frag_frag,Q_Electrostatic
|
||||
|
||||
|
||||
E = 0.0d0
|
||||
|
||||
do w1=1,nfrag-1
|
||||
ks1 = indx_frag_start(w1)
|
||||
n1 = size_frag(w1)
|
||||
Rin1 = switch_Rin(kfrag(w1))
|
||||
Rout1 = switch_Rout(kfrag(w1))
|
||||
call pspw_Q_cm(n1,ks1,amass,rion,rw1_cm)
|
||||
|
||||
do w2=w1+1,nfrag
|
||||
ks2 = indx_frag_start(w2)
|
||||
n2 = size_frag(w2)
|
||||
Rin2 = switch_Rin(kfrag(w2))
|
||||
Rout2 = switch_Rout(kfrag(w2))
|
||||
call pspw_Q_cm(n2,ks2,amass,rion,rw2_cm)
|
||||
|
||||
Rin = 0.5d0*(Rin1 +Rin2)
|
||||
Rout = 0.5d0*(Rout1+Rout2)
|
||||
E = E + pspw_Q_E_frag_frag(n1,ks1,rw1_cm,
|
||||
> n2,ks2,rw2_cm,
|
||||
> Rin,Rout,
|
||||
> katm,psp_zv,rion,
|
||||
> nshl3d,rcell)
|
||||
end do
|
||||
end do
|
||||
|
||||
c **** take out MM/MM Coulomb self energy ****
|
||||
do w1 = 1,nfrag
|
||||
ks1 = indx_frag_start(w1)
|
||||
do a=1,size_frag(w1)
|
||||
do b=a+1,size_frag(w1)
|
||||
kk1 = ks1+a-1
|
||||
kk2 = ks1+b-1
|
||||
E = E - Q_Electrostatic(rion(1,kk1),psp_zv(katm(kk1)),
|
||||
> rion(1,kk2),psp_zv(katm(kk2)))
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
pspw_Q_E = E
|
||||
return
|
||||
end
|
||||
|
||||
c *********************************************
|
||||
c * *
|
||||
c * pspw_Q_E_frag_frag *
|
||||
c * *
|
||||
c *********************************************
|
||||
|
||||
real*8 function pspw_Q_E_frag_frag(n1,ks1,rw1_cm,
|
||||
> n2,ks2,rw2_cm,
|
||||
> Rin,Rout,
|
||||
> katm,psp_zv,rion,
|
||||
> nshl3d,rcell)
|
||||
implicit none
|
||||
integer n1,ks1
|
||||
real*8 rw1_cm(3)
|
||||
integer n2,ks2
|
||||
real*8 rw2_cm(3)
|
||||
real*8 Rin,Rout
|
||||
integer katm(*)
|
||||
real*8 psp_zv(*)
|
||||
real*8 rion(3,*)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
|
||||
|
||||
* **** local variables ****
|
||||
integer l,a,b,kk1,kk2
|
||||
real*8 x,y,z,rcm,Sm
|
||||
real*8 r1(3),q1
|
||||
real*8 r2(3),q2
|
||||
real*8 E
|
||||
|
||||
* **** external functions ****
|
||||
real*8 pspw_Q_Switching,Q_Electrostatic
|
||||
external pspw_Q_Switching,Q_Electrostatic
|
||||
|
||||
E = 0.0d0
|
||||
do l=1,nshl3d
|
||||
x = rw1_cm(1) - rw2_cm(1) + rcell(l,1)
|
||||
y = rw1_cm(2) - rw2_cm(2) + rcell(l,2)
|
||||
z = rw1_cm(3) - rw2_cm(3) + rcell(l,3)
|
||||
rcm = dsqrt(x**2 + y**2 + z**2)
|
||||
if (rcm.lt.Rout) then
|
||||
Sm = pspw_Q_Switching(Rin,Rout,rcm)-1.0d0
|
||||
kk1 = ks1
|
||||
do a = 1,n1
|
||||
r1(1) = rion(1,kk1) + rcell(l,1)
|
||||
r1(2) = rion(2,kk1) + rcell(l,2)
|
||||
r1(3) = rion(3,kk1) + rcell(l,3)
|
||||
q1 = psp_zv(katm(kk1))
|
||||
|
||||
kk2 = ks2
|
||||
do b = 1,n2
|
||||
r2(1) = rion(1,kk2)
|
||||
r2(2) = rion(2,kk2)
|
||||
r2(3) = rion(3,kk2)
|
||||
q2 = psp_zv(katm(kk2))
|
||||
E = E + Sm*Q_Electrostatic(r1,q1,r2,q2)
|
||||
kk2 = kk2 + 1
|
||||
end do
|
||||
kk1 = kk1 + 1
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
|
||||
pspw_Q_E_frag_frag = E
|
||||
return
|
||||
end
|
||||
|
||||
* ******************************
|
||||
* * *
|
||||
* * pspw_Q_Switching *
|
||||
* * *
|
||||
* ******************************
|
||||
real*8 function pspw_Q_Switching(Rin,Rout,r)
|
||||
implicit none
|
||||
real*8 Rin,Rout,r
|
||||
|
||||
* **** ST2 parameters ****
|
||||
c real*8 Rl,Ru,c1,c2
|
||||
c parameter (Rl=2.0160d0/0.529177d0)
|
||||
c parameter (Ru=3.1287d0/0.529177d0)
|
||||
c parameter (c1 = (Ru-Rl)**3)
|
||||
c parameter (c2 = (3.0d0*Ru-Rl))
|
||||
real*8 c1,c2,s
|
||||
|
||||
if (r.le.Rin) then
|
||||
s = 0.0d0
|
||||
else if (r.ge.Rout) then
|
||||
s = 1.0d0
|
||||
else
|
||||
c1 = (Rout-Rin)**3
|
||||
c2 = 3.0d0*Rout-Rin
|
||||
s = (r-Rin)**2 * (c2-2.0d0*r)/c1
|
||||
end if
|
||||
|
||||
pspw_Q_Switching = s
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ******************************
|
||||
* * *
|
||||
* * pspw_Q_cm *
|
||||
* * *
|
||||
* ******************************
|
||||
subroutine pspw_Q_cm(n,ks,amass,rion,rcm)
|
||||
implicit none
|
||||
integer n,ks
|
||||
real*8 amass(*),rion(3,*),rcm(3)
|
||||
|
||||
* **** local variables ****
|
||||
integer k,kk
|
||||
real*8 m
|
||||
|
||||
m = 0.0d0
|
||||
rcm(1) = 0.0d0
|
||||
rcm(2) = 0.0d0
|
||||
rcm(3) = 0.0d0
|
||||
kk = ks
|
||||
do k=1,n
|
||||
rcm(1) = rcm(1) + amass(kk)*rion(1,kk)
|
||||
rcm(2) = rcm(2) + amass(kk)*rion(2,kk)
|
||||
rcm(3) = rcm(3) + amass(kk)*rion(3,kk)
|
||||
m = m + amass(kk)
|
||||
kk = kk + 1
|
||||
end do
|
||||
rcm(1) = rcm(1)/m
|
||||
rcm(2) = rcm(2)/m
|
||||
rcm(3) = rcm(3)/m
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
c
|
||||
c $Id: pspw_Q.F,v 1.12 2007-03-22 20:46:23 bylaska Exp $
|
||||
c
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
* **********************************
|
||||
* * *
|
||||
* * Q_Electrostatic_mForce *
|
||||
* * *
|
||||
* **********************************
|
||||
|
||||
* This routine calculates minus the Electrostatic force between two centers
|
||||
*
|
||||
subroutine Q_Electrostatic_mForce(r1,q1,f1,r2,q2,f2)
|
||||
implicit none
|
||||
real*8 r1(3),q1
|
||||
real*8 f1(3)
|
||||
real*8 r2(3),q2
|
||||
real*8 f2(3)
|
||||
|
||||
* **** local variables ****
|
||||
real*8 xo,yo,zo,d0,d2
|
||||
real*8 der
|
||||
real*8 fx,fy,fz
|
||||
|
||||
|
||||
xo = r2(1)-r1(1)
|
||||
yo = r2(2)-r1(2)
|
||||
zo = r2(3)-r1(3)
|
||||
|
||||
d2 = (xo*xo + yo*yo + zo*zo)
|
||||
d0 = dsqrt(d2)
|
||||
|
||||
der = -q1*q2/d2
|
||||
fx = -(xo/d0)*der
|
||||
fy = -(yo/d0)*der
|
||||
fz = -(zo/d0)*der
|
||||
|
||||
f2(1) = f2(1) - fx
|
||||
f2(2) = f2(2) - fy
|
||||
f2(3) = f2(3) - fz
|
||||
|
||||
f1(1) = f1(1) + fx
|
||||
f1(2) = f1(2) + fy
|
||||
f1(3) = f1(3) + fz
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
c *************************************
|
||||
c * *
|
||||
c * pspw_Q_fion *
|
||||
c * *
|
||||
c *************************************
|
||||
subroutine pspw_Q_fion(nion,nion_qm,katm,psp_zv,amass,
|
||||
> nfrag,indx_frag_start,size_frag,
|
||||
> kfrag,
|
||||
> switch_Rin,switch_Rout,
|
||||
> nshl3d,rcell,
|
||||
> rion,fion)
|
||||
implicit none
|
||||
integer nion,nion_qm
|
||||
integer katm(*)
|
||||
real*8 psp_zv(*),amass(*)
|
||||
integer nfrag
|
||||
integer indx_frag_start(*),size_frag(*)
|
||||
integer kfrag(*)
|
||||
real*8 switch_Rin(*),switch_Rout(*)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
real*8 rion(3,*),fion(3,*)
|
||||
|
||||
|
||||
* **** local variables ****
|
||||
integer w1,a,ks1,kk1,n1
|
||||
integer w2,b,ks2,kk2,n2
|
||||
real*8 Rin1,Rout1,rw1_cm(3),m1
|
||||
real*8 Rin2,Rout2,rw2_cm(3),m2
|
||||
real*8 Rin,Rout
|
||||
|
||||
* **** external functions ****
|
||||
real*8 pspw_Q_E_frag_frag
|
||||
external pspw_Q_E_frag_frag
|
||||
|
||||
|
||||
|
||||
do w1=1,nfrag-1
|
||||
ks1 = indx_frag_start(w1)
|
||||
n1 = size_frag(w1)
|
||||
Rin1 = switch_Rin(kfrag(w1))
|
||||
Rout1 = switch_Rout(kfrag(w1))
|
||||
call pspw_Q_cmm(n1,ks1,amass,rion,rw1_cm,m1)
|
||||
|
||||
do w2=w1+1,nfrag
|
||||
ks2 = indx_frag_start(w2)
|
||||
n2 = size_frag(w2)
|
||||
Rin2 = switch_Rin(kfrag(w2))
|
||||
Rout2 = switch_Rout(kfrag(w2))
|
||||
call pspw_Q_cmm(n2,ks2,amass,rion,rw2_cm,m2)
|
||||
|
||||
Rin = 0.5d0*(Rin1 +Rin2)
|
||||
Rout = 0.5d0*(Rout1+Rout2)
|
||||
call pspw_Q_fion_frag_frag(n1,ks1,rw1_cm,m1,
|
||||
> n2,ks2,rw2_cm,m2,
|
||||
> Rin,Rout,
|
||||
> katm,psp_zv,amass,rion,fion,
|
||||
> nshl3d,rcell)
|
||||
end do
|
||||
end do
|
||||
|
||||
c **** take out MM/MM Coulomb self energy ****
|
||||
do w1 = 1,nfrag
|
||||
ks1 = indx_frag_start(w1)
|
||||
do a=1,size_frag(w1)
|
||||
do b=a+1,size_frag(w1)
|
||||
kk1 = ks1+a-1
|
||||
kk2 = ks1+b-1
|
||||
call Q_Electrostatic_mForce(
|
||||
> rion(1,kk1),psp_zv(katm(kk1)),fion(1,kk1),
|
||||
> rion(1,kk2),psp_zv(katm(kk2)),fion(1,kk2))
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
c *********************************************
|
||||
c * *
|
||||
c * pspw_Q_fion_frag_frag *
|
||||
c * *
|
||||
c *********************************************
|
||||
|
||||
subroutine pspw_Q_fion_frag_frag(n1,ks1,rw1_cm,m1,
|
||||
> n2,ks2,rw2_cm,m2,
|
||||
> Rin,Rout,
|
||||
> katm,psp_zv,amass,rion,fion,
|
||||
> nshl3d,rcell)
|
||||
implicit none
|
||||
integer n1,ks1
|
||||
real*8 rw1_cm(3),m1
|
||||
integer n2,ks2
|
||||
real*8 rw2_cm(3),m2
|
||||
real*8 Rin,Rout
|
||||
integer katm(*)
|
||||
real*8 psp_zv(*)
|
||||
real*8 amass(*)
|
||||
real*8 rion(3,*)
|
||||
real*8 fion(3,*)
|
||||
integer nshl3d
|
||||
real*8 rcell(nshl3d,3)
|
||||
|
||||
|
||||
* **** local variables ****
|
||||
integer l,a,b,kk1,kk2
|
||||
real*8 x,y,z,rcm,Sm,dSm
|
||||
real*8 r1(3),q1,f1(3)
|
||||
real*8 r2(3),q2,f2(3)
|
||||
real*8 E
|
||||
|
||||
* **** external functions ****
|
||||
real*8 Q_Electrostatic
|
||||
external Q_Electrostatic
|
||||
|
||||
E = 0.0d0
|
||||
do l=1,nshl3d
|
||||
x = rw1_cm(1) - rw2_cm(1) + rcell(l,1)
|
||||
y = rw1_cm(2) - rw2_cm(2) + rcell(l,2)
|
||||
z = rw1_cm(3) - rw2_cm(3) + rcell(l,3)
|
||||
rcm = dsqrt(x**2 + y**2 + z**2)
|
||||
if (rcm.lt.Rout) then
|
||||
call pspw_Q_dSwitching(Rin,Rout,rcm,Sm,dSm)
|
||||
Sm = Sm - 1.0d0
|
||||
|
||||
* *** calculate E, and -Sm*grad(E) ****
|
||||
E = 0.0d0
|
||||
kk1 = ks1
|
||||
do a = 1,n1
|
||||
r1(1) = rion(1,kk1) + rcell(l,1)
|
||||
r1(2) = rion(2,kk1) + rcell(l,2)
|
||||
r1(3) = rion(3,kk1) + rcell(l,3)
|
||||
q1 = psp_zv(katm(kk1))
|
||||
|
||||
kk2 = ks2
|
||||
do b = 1,n2
|
||||
r2(1) = rion(1,kk2)
|
||||
r2(2) = rion(2,kk2)
|
||||
r2(3) = rion(3,kk2)
|
||||
q2 = psp_zv(katm(kk2))
|
||||
E = E + Sm*Q_Electrostatic(r1,q1,r2,q2)
|
||||
|
||||
call Q_Electrostatic_Force(r1,q1,f1,r2,q2,f2)
|
||||
fion(1,kk1) = fion(1,kk1) + Sm*f1(1)
|
||||
fion(2,kk1) = fion(2,kk1) + Sm*f1(2)
|
||||
fion(3,kk1) = fion(3,kk1) + Sm*f1(3)
|
||||
fion(1,kk2) = fion(1,kk2) + Sm*f2(1)
|
||||
fion(2,kk2) = fion(2,kk2) + Sm*f2(2)
|
||||
fion(3,kk2) = fion(3,kk2) + Sm*f2(3)
|
||||
|
||||
kk2 = kk2 + 1
|
||||
end do
|
||||
kk1 = kk1 + 1
|
||||
end do
|
||||
|
||||
* *** calculate -E*grad(S) ****
|
||||
kk1 = ks1
|
||||
do a=1,n1
|
||||
fion(1,kk1) = fion(1,kk1) - E*dSm*(x/rcm)*amass(kk1)/m1
|
||||
fion(2,kk1) = fion(1,kk1) - E*dSm*(y/rcm)*amass(kk1)/m1
|
||||
fion(3,kk1) = fion(1,kk1) - E*dSm*(z/rcm)*amass(kk1)/m1
|
||||
kk1 = kk1 + 1
|
||||
end do
|
||||
kk2 = ks2
|
||||
do b=1,n2
|
||||
fion(1,kk2) = fion(1,kk2) + E*dSm*(x/rcm)*amass(kk2)/m2
|
||||
fion(2,kk2) = fion(1,kk2) + E*dSm*(y/rcm)*amass(kk2)/m2
|
||||
fion(3,kk2) = fion(1,kk2) + E*dSm*(z/rcm)*amass(kk2)/m2
|
||||
kk2 = kk2 + 1
|
||||
end do
|
||||
|
||||
end if
|
||||
end do
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
* ******************************
|
||||
* * *
|
||||
* * pspw_Q_dSwitching *
|
||||
* * *
|
||||
* ******************************
|
||||
subroutine pspw_Q_dSwitching(Rin,Rout,r,S,dS)
|
||||
implicit none
|
||||
real*8 Rin,Rout,r,S,dS
|
||||
|
||||
* **** ST2 parameters ****
|
||||
c real*8 Rl,Ru,c1,c2
|
||||
c parameter (Rl=2.0160d0/0.529177d0)
|
||||
c parameter (Ru=3.1287d0/0.529177d0)
|
||||
c parameter (c1 = (Ru-Rl)**3)
|
||||
c parameter (c2 = (3.0d0*Ru-Rl))
|
||||
real*8 c1,c2
|
||||
|
||||
if (r.le.Rin) then
|
||||
S = 0.0d0
|
||||
dS = 0.0d0
|
||||
else if (r.ge.Rout) then
|
||||
S = 1.0d0
|
||||
dS = 0.0d0
|
||||
else
|
||||
c1 = (Rout-Rin)**3
|
||||
c2 = 3.0d0*Rout-Rin
|
||||
S = (r-Rin)**2 * (c2-2.0d0*r)/c1
|
||||
dS = 2.0d0*(r-Rin) * (c2-2.0d0*r)/c1
|
||||
> - 2.0d0*(r-Rin)**2 /c1
|
||||
|
||||
end if
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
* ******************************
|
||||
* * *
|
||||
* * pspw_Q_cmm *
|
||||
* * *
|
||||
* ******************************
|
||||
subroutine pspw_Q_cmm(n,ks,amass,rion,rcm,m)
|
||||
implicit none
|
||||
integer n,ks
|
||||
real*8 amass(*),rion(3,*),rcm(3)
|
||||
real*8 m
|
||||
|
||||
* **** local variables ****
|
||||
integer k,kk
|
||||
|
||||
m = 0.0d0
|
||||
rcm(1) = 0.0d0
|
||||
rcm(2) = 0.0d0
|
||||
rcm(3) = 0.0d0
|
||||
kk = ks
|
||||
do k=1,n
|
||||
rcm(1) = rcm(1) + amass(kk)*rion(1,kk)
|
||||
rcm(2) = rcm(2) + amass(kk)*rion(2,kk)
|
||||
rcm(3) = rcm(3) + amass(kk)*rion(3,kk)
|
||||
m = m + amass(kk)
|
||||
kk = kk + 1
|
||||
end do
|
||||
rcm(1) = rcm(1)/m
|
||||
rcm(2) = rcm(2)/m
|
||||
rcm(3) = rcm(3)/m
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load diff
|
|
@ -1,31 +1,18 @@
|
|||
c
|
||||
c $Id: pspw_qmmm.fh,v 1.6 2002-02-26 22:45:56 bylaska Exp $
|
||||
c $Id: pspw_qmmm.fh,v 1.7 2007-03-22 20:46:22 bylaska Exp $
|
||||
c
|
||||
|
||||
* ***** pspw_qmmm common block ************************************
|
||||
real*8 ekw,ekw_total
|
||||
integer r2(2),r1(2),r0(2),mass(2)
|
||||
integer katm(2)
|
||||
integer natm(2)
|
||||
integer atom(2)
|
||||
integer frag_from_indx(2)
|
||||
integer indx_from_frag(2),dsqm_frag(2)
|
||||
integer size_frag(2),shake_frag(2)
|
||||
integer ekw_count
|
||||
integer nion,nkatm,nfrag,max_frag_length
|
||||
integer qmmm_type,geom
|
||||
common / pspw_qmmm / ekw,ekw_total,r2,r1,r0,mass,
|
||||
> katm, ! katm list
|
||||
> natm, ! natm list
|
||||
> atom, ! atom list
|
||||
> frag_from_indx,
|
||||
> indx_from_frag,dsqm_frag,
|
||||
> size_frag,shake_frag,
|
||||
> ekw_count,
|
||||
> nion, ! number QM/MM ions
|
||||
> nkatm, ! number of kind of QM/MM ions
|
||||
> nfrag, ! number of fragments
|
||||
> max_frag_length, ! maximum fragment size
|
||||
> qmmm_type, ! 1-model potential,2-switching charge,3-model charge
|
||||
> geom
|
||||
integer indx_frag_start(2)
|
||||
integer size_frag(2)
|
||||
integer kfrag(2)
|
||||
integer switch_Rin(2),switch_Rout(2)
|
||||
integer nshl3d,rcell(2)
|
||||
integer nfrag,nkfrag
|
||||
logical qmmm_found
|
||||
common / pspw_qmmm / indx_frag_start,
|
||||
> size_frag,
|
||||
> kfrag,switch_Rin,switch_Rout,
|
||||
> rcell,nshl3d,
|
||||
> nfrag,nkfrag,qmmm_found
|
||||
* ******************************************************************
|
||||
|
|
|
|||
|
|
@ -9,26 +9,40 @@
|
|||
#include "util.fh"
|
||||
c
|
||||
c
|
||||
integer ind ! Index of matched directive
|
||||
integer num_dirs ! No. of known directives
|
||||
parameter (num_dirs = 6)
|
||||
|
||||
character*30 dirs(num_dirs)
|
||||
character*255 test
|
||||
|
||||
integer num_dirs ! No. of known directives
|
||||
parameter (num_dirs = 7)
|
||||
character*30 dirs(num_dirs)
|
||||
data dirs / 'lj_ion_parameters:',
|
||||
> 'lj_mm_parameters:',
|
||||
> 'polarization_mm_parameter:',
|
||||
> 'sr_basis:',
|
||||
> 'lr_cutoff:',
|
||||
> 'fragment',
|
||||
> 'end'/
|
||||
|
||||
integer num_fdirs ! No. of known directives
|
||||
parameter (num_fdirs = 3)
|
||||
character*30 fdirs(num_fdirs)
|
||||
data fdirs / 'size',
|
||||
> 'index_start',
|
||||
> 'end'/
|
||||
c
|
||||
logical value
|
||||
character*4 element
|
||||
character*80 rtdb_name
|
||||
integer element_length,basis_length,tmp(2)
|
||||
integer ni,nkf,nfrag,frag_size,index_start(30)
|
||||
real*8 se(2),alpha
|
||||
character*255 test
|
||||
integer ind ! Index of matched directive
|
||||
|
||||
* **** external functions ****
|
||||
character*7 c_index_name
|
||||
external c_index_name
|
||||
|
||||
nkf = 0
|
||||
nfrag = 0
|
||||
|
||||
10 if (.not. inp_read())
|
||||
> call errquit(
|
||||
|
|
@ -41,7 +55,7 @@ c
|
|||
> 'pspw_qmmm_input: unknown directive', 0, INPUT_ERR)
|
||||
|
||||
|
||||
goto ( 100, 200, 300, 400, 500,
|
||||
goto ( 100, 200, 300, 400, 500, 600,
|
||||
> 9999) ind
|
||||
call errquit(
|
||||
> 'pspw_qmmm_input: unimplemented directive', ind, INPUT_ERR)
|
||||
|
|
@ -165,6 +179,71 @@ c
|
|||
|
||||
goto 10
|
||||
|
||||
c
|
||||
c fragment
|
||||
c
|
||||
600 nkf = nkf + 1
|
||||
rtdb_name = 'pspw_qmmm_nkfrag'
|
||||
value = rtdb_put(rtdb,rtdb_name,mt_int,1,nkf)
|
||||
if (.not. value)
|
||||
> call errquit(
|
||||
> 'pspw_qmmm_input:failed writing nkfrag',0, RTDB_ERR)
|
||||
|
||||
601 if (.not. inp_read())
|
||||
> call errquit('pspw_qmmm_input: inp_read failed',0,INPUT_ERR)
|
||||
if (.not. inp_a(test))
|
||||
>call errquit('pspw_qmmm_input:failed reading keyword',0,INPUT_ERR)
|
||||
if (.not. inp_match(num_fdirs, .false.,test,fdirs,ind))
|
||||
> call errquit('pspw_qmmm_input: unknown directive', 0, INPUT_ERR)
|
||||
|
||||
|
||||
goto ( 610, 620,
|
||||
> 699) ind
|
||||
call errquit('pspw_qmmm_input:directive not found',ind,INPUT_ERR)
|
||||
|
||||
c
|
||||
c size
|
||||
c
|
||||
610 if (.not.inp_i(frag_size)) frag_size = 1
|
||||
rtdb_name = 'pspw_qmmm_frag_size:'//c_index_name(nkf)
|
||||
value = rtdb_put(rtdb,rtdb_name,mt_int,1,frag_size)
|
||||
if (.not. value)
|
||||
> call errquit(
|
||||
> 'pspw_qmmm_input:failed writing frag_size',0, RTDB_ERR)
|
||||
|
||||
goto 601
|
||||
c
|
||||
c index_start
|
||||
c
|
||||
620 ni = 0
|
||||
do while (inp_i(index_start(ni+1)))
|
||||
ni = ni + 1
|
||||
end do
|
||||
rtdb_name = 'pspw_qmmm_frag_index_start:'//c_index_name(nkf)
|
||||
value = rtdb_put(rtdb,rtdb_name,mt_int,ni,index_start)
|
||||
if (.not. value)
|
||||
> call errquit(
|
||||
> 'pspw_qmmm_input:failed writing index_start',0, RTDB_ERR)
|
||||
|
||||
rtdb_name = 'pspw_qmmm_frag_nindex_start:'//c_index_name(nkf)
|
||||
value = rtdb_put(rtdb,rtdb_name,mt_int,1,ni)
|
||||
if (.not. value)
|
||||
> call errquit(
|
||||
> 'pspw_qmmm_input:failed writing nindex_start',0, RTDB_ERR)
|
||||
|
||||
nfrag = nfrag + ni
|
||||
rtdb_name = 'pspw_qmmm_nfrag'
|
||||
value = rtdb_put(rtdb,rtdb_name,mt_int,1,nfrag)
|
||||
if (.not. value)
|
||||
> call errquit(
|
||||
> 'pspw_qmmm_input:failed writing nfrag',0, RTDB_ERR)
|
||||
|
||||
goto 601
|
||||
|
||||
|
||||
699 continue
|
||||
goto 10
|
||||
|
||||
9999 continue
|
||||
|
||||
return
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue