1_polar_dipole.inp

* analysis for molec polarizabilities and dipole moments
*

ioformat extended
! Note: You need to change below QM dipole and alpha and put your QM results
! This script will read the resi_mp2.crd


set residue clet  
! set residue chlb

set TOPPAR .
stream toppar_drude_master_example.str
stream toppar_drude_example.str

! qm optimized @resi_mp.crd
set mp2dir ./gauss_opt

!non-bond parameters
set 3 999.0  ! cutim
set 4 999.0  ! cutnb
set 5 980.0  ! ctonnb
set 6 990.0  ! ctofnb
set 7 switch
set 8 atom
set 9 vatom

read sequence @residue 1
generate @residue first none last none setup warn drude dmass 0.4
bomlev 0

open read card unit 21 name @mp2dir/@residue_mp2.crd
read coor card unit 21

coor sdrude
coor shake
coor copy comp

update cdie inbfrq -1 ihbfrq 0 -
@7 @8 @9 vswitch cutnb @4 ctofnb @6 ctonnb @5

cons harm force 10000000.0 sele .not. type D* show end
cons harm force 10000000.0 sele type dum show end
MINI SD nstep 200 tolgrd 0.00001
MINI ABNR nstep 500 tolgrd 0.00001
cons harm force 0.0 sele all end

MINI SD nstep 200 tolgrd 0.00001
MINI ABNR nstep 500 tolgrd 0.00001

! start analysis
open unit 50 write form name polar.dipole.@residue.dat
coor dipole oxyz select segid @residue end
echu 50
echo QM dipole:
echo mp2/cc-pvQZ dipole: X= -1.1386 Y=1.7948 Z=0.0000  Tot=2.1254 
echo mm dipole: X=  ?xdip  Y=  ?ydip  Z=  ?zdip  Tot= ?rdip
echo 

stream polar.efield.str
echo QM alpha:
echo mp2/cc-pvQZ alpha: XX= 6.249 YY= 6.359 ZZ= 4.962 Tot= 5.857
echo mm alpha: xx= @alphaxx yy= @alphayy  zz=  @alphazz  tot= @malpha

stop