-
Notifications
You must be signed in to change notification settings - Fork 19
scripting
import os
import numpy as np
import matplotlib.pyplot as plt
import elements
import basUtils
def makeclean( ):
import os
[ os.remove(f) for f in os.listdir(".") if f.endswith(".so") ]
[ os.remove(f) for f in os.listdir(".") if f.endswith(".o") ]
[ os.remove(f) for f in os.listdir(".") if f.endswith(".pyc") ]
makeclean( ) # force to recompile
import ProbeParticle as PP
PP.loadParams( 'params.ini' ) # load parametes from ini file
atoms = basUtils.loadAtoms('input.xyz', elements.ELEMENT_DICT )
Rs = np.array([atoms[1],atoms[2],atoms[3]]);
iZs = np.array( atoms[0])
if not PP.params['PBC' ]:
PP.autoGeom( Rs, shiftXY=True, fitCell=True, border=3.0 )
Rs[0] += PP.params['moleculeShift' ][0] # shift molecule so that we sample reasonable part of potential
Rs[1] += PP.params['moleculeShift' ][1]
Rs[2] += PP.params['moleculeShift' ][2]
Rs = np.transpose( Rs, (1,0) ).copy()
Qs = np.array( atoms[4] )
if PP.params['PBC' ]:
iZs,Rs,Qs = PP.PBCAtoms( iZs, Rs, Qs, avec=PP.params['gridA'], bvec=PP.params['gridB'] )
### Define and allocate arrays
do this before simulation, in case it will crash
dz = PP.params['scanStep'][2] zTips = np.arange( PP.params['scanMin'][2], PP.params['scanMax'][2]+0.00001, dz )[::-1]; ntips = len(zTips); print " zTips : ",zTips rTips = np.zeros((ntips,3)) rs = np.zeros((ntips,3)) fs = np.zeros((ntips,3))
rTips[:,0] = 1.0 rTips[:,1] = 1.0 rTips[:,2] = zTips
PP.setTip()
xTips = np.arange( PP.params['scanMin'][0], PP.params['scanMax'][0]+0.00001, 0.1 ) yTips = np.arange( PP.params['scanMin'][1], PP.params['scanMax'][1]+0.00001, 0.1 ) extent=( xTips[0], xTips[-1], yTips[0], yTips[-1] ) fzs = np.zeros(( len(zTips), len(yTips ), len(xTips ) ));
nslice = 10;
FFparams = PP.loadSpecies ( 'atomtypes.ini' ) C6,C12 = PP.getAtomsLJ( PP.params['probeType'], iZs, FFparams )
print " # ============ define Grid "
cell =np.array([ PP.params['gridA'], PP.params['gridB'], PP.params['gridC'], ]).copy()
gridN = PP.params['gridN']
FF = np.zeros( (gridN[2],gridN[1],gridN[0],3) )