Geometry optimization of NaCl clusters
Use this short script to drive CP2K
#!/bin/bash --login
#PBS -N cp2k
#PBS -l select=1
#PBS -l walltime=0:20:0
#PBS -A y14
cd $PBS_O_WORKDIR
module load cp2k
for ii in 2 4 6 8 10 12
do
sed -e "s/MY_SUPERCELL/${ii}/g" template.inp > input_${ii}.inp
aprun -n 2 cp2k.popt input_${ii}.inp > NaCl_supercell_${ii}.out
done
where the template input template.inp is this geometry optimization using the classical forcefield FIST module of CP2K.
@SET NREP MY_SUPERCELL
@SET OPTIMIZER LBFGS # BFGS
&FORCE_EVAL
METHOD Fist
&MM
&FORCEFIELD
&CHARGE
ATOM Na
CHARGE +1.000
&END CHARGE
&CHARGE
ATOM Cl
CHARGE -1.000
&END CHARGE
&NONBONDED
&BMHFT
map_atoms NA NA
atoms NA NA
RCUT 10.0
&END BMHFT
&BMHFT
map_atoms NA CL
atoms NA CL
RCUT 10.0
&END BMHFT
&BMHFT
map_atoms CL CL
atoms CL CL
RCUT 10.0
&END BMHFT
&END NONBONDED
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE spme
ALPHA .35
GMAX 12*${NREP}
O_SPLINE 6
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
#ABC 5.620 5.620 5.620
ABC 2*5.620 2*5.620 2*5.620
MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
&END CELL
&TOPOLOGY
COORD_FILE_NAME NaCl.pdb
COORDINATE PDB
CONN_FILE_FORMAT OFF
MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT NaCl
RUN_TYPE GEO_OPT
&END GLOBAL
&MOTION
&GEO_OPT
OPTIMIZER ${OPTIMIZER}
&END
&END MOTION
You will also need to initial geometry file NaCl.pdb:
REMARK 4 NaCl COMPLIES WITH FORMAT V. 2.0 CRYST1 5.620 5.620 5.620 90.00 90.00 90.00 P 1 HETATM 1 Na NMO 0 0.000 0.000 0.000 1.00 0.00 NMO Na HETATM 2 Cl CMO 0 2.810 2.810 2.810 1.00 0.00 CMO Cl HETATM 3 Na NMO 0 0.000 2.810 2.810 1.00 0.00 NMO Na HETATM 4 Cl CMO 0 2.810 0.000 0.000 1.00 0.00 CMO Cl HETATM 5 Na NMO 0 2.810 0.000 2.810 1.00 0.00 NMO Na HETATM 6 Cl CMO 0 0.000 2.810 0.000 1.00 0.00 CMO Cl HETATM 7 Na NMO 0 2.810 2.810 0.000 1.00 0.00 NMO Na HETATM 8 Cl CMO 0 0.000 0.000 2.810 1.00 0.00 CMO Cl TER 9 UNK 0 END
The script runs very quickly when the LBFGS optimizer is used. But see what happens if we switch to the BFGS optimizer instead (change the OPTIMIZER variable in the cp2k input file - you might want to reduce the size of the supercells in the driver script - NREP varying from 1 to 6 perhaps). Look at the timings that cp2k prints at the end of a run and see if you can see the culprit. Also look for warnings in your outputs (this is a good habit to get into).
How does the conjugate gradients optimizer compare to LBFGS in efficiency for this system?
Cell optimization of NaCl
For studying many properties of solid materials it is important that the lattice parameters used in a simulation are close to equilibrium for the model chemistry (Hamiltonian) used. Otherwise, large stresses can be present that complicate comparison to experiment. Successful cell optimization requires that the energy changes smoothly with cell volume - and for this the energy cutoff is the most important parameter. The input file template below template.inp can be used with driver script to examine how the energy volume curve of NaCl changes with the PW cutoff.
@SET SCALE_FACTOR MY_SCALING
@SET NREP 1
@SET OPTIMIZER BFGS # LBFGS
@SET CUTOFF MY_CUTOFF
@SET SAFTEY_CUTOFF 1.1
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF ${CUTOFF}
REL_CUTOFF 60
&END MGRID
&QS
EPS_DEFAULT 1.0E-12
&END QS
&SCF
SCF_GUESS RESTART
&OT ON
MINIMIZER DIIS
&END OT
&END SCF
&XC
&XC_FUNCTIONAL Pade
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 5.620*${SCALE_FACTOR} 5.620*${SCALE_FACTOR} 5.620*${SCALE_FACTOR}
MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
#&CELL_REF
# ABC 5.620*${SAFETY_FACTOR} 5.620*${SAFETY_FACTOR} 5.620*${SAFETY_FACTOR}
# MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
#&END
&END CELL
&COORD
scaled
Na 0.000 0.000 0.000
Cl 0.500 0.500 0.500
Na 0.000 0.500 0.500
Cl 0.500 0.000 0.000
Na 0.500 0.000 0.500
Cl 0.000 0.500 0.000
Na 0.500 0.500 0.000
Cl 0.000 0.000 0.500
&END
&TOPOLOGY
MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
&END TOPOLOGY
&KIND Na
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PADE-q9
&END KIND
&KIND Cl
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PADE-q7
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT NaCl
RUN_TYPE ENERGY
&END GLOBAL
&MOTION
&GEO_OPT
OPTIMIZER ${OPTIMIZER}
&END
&CELL_OPT
KEEP_SYMMETRY
&END
&END MOTION
The driver script could be something like
#!/bin/bash --login
#PBS -N cp2k
#PBS -l select=1
#PBS -l walltime=0:20:0
#PBS -A y14
cd $PBS_O_WORKDIR
module load cp2k
CUTOFF="280"
for ii in 0.800 0.825 0.850 0.875 0.900 0.925 0.950 0.975 1.000 1.025 1.050 1.075 1.100
do
#this assumes that the input template is in NaCl_QS.inp
sed -e "s/MY_SCALING/${ii}/g" NaCl_QS.inp > temp.inp
sed -e "s/MY_CUTOFF/${CUTOFF}/g" temp.inp > input_${ii}.inp
#this line should changed to point to your cp2k executable
aprun -n 2 cp2k.popt input_${ii}.inp > NaCl_${CUTOFF}_${ii}.out
done
You can extract energies from the outputs with a command like
grep 'ENERGY|' *out | awk '{print $10}' > NaCl_energy_volume.dat
and you can plot the results in your favourite graphing software.
What is happening here? Try changing the PW cutoff (defined in the driver script) and using the CELL_REF variable.
Copy the input template to a new file and change the RUN_TYPE to CELL_OPT. You'll also need to ask the code to calculate the stress tensor (in FORCE_EVAL section) ANALYTICALly! Also define the CUTOFF and SCALING_FACTOR. Start the cell optimization from small cells (SCALING_FACTOR 0.85) or large cells (SCALING_FACTOR 1.10) - do you get the same results?
If you are running on ARCHER, use some larger parallel jobs to check whether increasing the supercell size (NREP variable) affects the results.