Scan job error at cartesian coordinates

2 posts / 0 new
Last post
xixhaha
Scan job error at cartesian coordinates

Hi everyone:
When doing the scan calculation of QM/MM calculation, I usually meet a problem about "Conversion of residue 1 failed , HDLC failure gauge: 2000".
The residue 1 is in my QM region, and I check the geometry but can not see anything strange.
I read the relevant questions that Tom answered on the forum. He said that to solve this kind of problem, we should use Cartesian coordinates and froze the atoms at the same time.
I try to restrict the atom 266 and atom 4402 (two atoms in the scan process), but the two atoms can not be restrainted and back to the reactant species.
I think there is something wrong with the input script, I really need everyone's help. Can everyone help me check and change my input script?

==================== =================== ==================== ======================================== ===================
MM Energies (kcal/mol) total : -7.40357489E+04
short range: 6.59150532E+03 coulombic : -8.91566442E+04 3body : 0.00000000E+00
bond : 4.74680054E+03 angle : 2.08153658E+03 str-bend : 0.00000000E+00
per dihed : 1.49220052E+03 harm dihed : 9.00608667E+01 angle-angle: 0.00000000E+00
ang-ang-tor: 0.00000000E+00 urey-brdly : 1.18791513E+02 inversion : 0.00000000E+00
Contribution to energy from turbomole: -16165.956374 (a.u.)
Contribution to energy from dl_poly: -117.980631 (a.u.)
Contribution to energy from additional MM energy terms: 0.000000 (a.u.)
-----------------------------------------------------------------------------------
QM/MM Energy: -16283.937005 (a.u.)
-----------------------------------------------------------------------------------
Energy calculation finished, energy: -1.628393701E+04

Converting Cartesians to HDLC
Testing convergence in cycle 2
Energy 2.1752E-06 Target: 1.0000E-05 converged? yes
Max step 6.0042E-03 Target: 3.6000E-03 converged? no component 40
RMS step 3.3605E-04 Target: 2.4000E-03 converged? yes
Max grad 5.4108E-03 Target: 9.0000E-04 converged? no component 40
RMS grad 3.0230E-04 Target: 6.0000E-04 converged? yes
Predicted step length 2.0000E-02
Trust radius 1.0000E-01

Converting HDLC to Cartesians
Conversion of residue 1 failed , HDLC failure gauge: 2000

Cyclic failure at residue 1, stopping

DL-FIND ERROR:
Residue conversion error

dlf_error called
jumped out of DL-Find with code 1
exit -1
exited DL-Find search: iret=-1
dl-find/hybrid.kill/============================= Tstep: 32.2 Ttot: 1050.4 ==
...hybrid.kill/dl_poly.kill/===================== Tstep: 0.0 Ttot: 1050.4 ==
MM termination code: 0
Energy and gradient evaluations:
QM energies 2
QM gradients 2
MM energies 2
MM gradients 2
0
ChemShell exiting code 0
the following objects are still declared:
<tag> <type> <links> <list> <name>
*hybrid.turbomole.coords fragment 1 main molecule 5
*dl-find.coo fragment 1 main CHARMM import
*dl-find.energy matrix 1 main Scratch energies for dl-find interface
*dl-find.gradient matrix 1 main Scratch gradient for dl-find interface
writing persistent objects to disk..
done
==================== =================== ==================== ======================================== ===================

my input script as below

==================== =================== ==================== ======================================== ===================
global sys_name_id
set sys_name_id sacc

set control_input_settings [ open control_input.${sys_name_id} w]
puts $control_input_settings " the jobname is ${sys_name_id}"
#
# define qmatoms
#
set qmatoms {258-268,453-467,4373-4429}

#define residues
#
set pdbresidues [ pdb_to_res "${sys_name_id}.pdb" ]
set myresidues [ inlist function=combine residues= $pdbresidues sets= {ENZY16 PROT302} target=fatone ]

#puts $myresidues

# define force field
#
set top { /home/wzlai/zcx/sacc/setup1/setup14/lib/top_all27_prot_lipid.rtf /home/wzlai/zcx/sacc/setup1/setup14/lib/top_sacc_prot.inp }
set prm /home/wzlai/zcx/sacc/setup1/setup14/lib/par_all27_prot_lipid_mod2e.prm

# load the connectivity:
#
load_connect_from_psf ${sys_name_id}.c ${sys_name_id}.psf

# read type, charges, group
#
source save_${sys_name_id}.chm

# read active atoms
#
set active_atom_numbers { 254-521 872-885 1089-1106 1160-1175 1245-1363 1549-1568 1724-1817 2067-2143 \
2986-3007 3024-3037 3428-3529 3876-3899 4050-4270 4293-4329 4373-4429 8198-8200 8603-8605 \
9026-9028 9665-9667 10799-10801 10865-10867 11144-11146 11159-11161 11225-11227 \
11246-11248 11552-11554 11648-11650 11807-11809 11816-11818 11930-11932 14468-14470 \
14495-14497 14948-14950 15113-15115 16310-16312 24716-24718 24866-24868 24914-24916 \
25193-25195 25208-25210 25211-25213 }

# define QM/MM settings:
#
set embedding_scheme shift
global qm_theory
set qm_theory turbomole
fragment hybrid.${qm_theory}.coords new persistent
set hybrid_args [ list \
coupling= $embedding_scheme \
cutoff= 99 \
groups = $groups \
atom_charges= $charges \
qm_region= $qmatoms \
qm_theory= $qm_theory : [ list \
read_control= yes \
scratchdir= /scr/wzlai \
scftype=uhf \
maxscratch= 4000 \
hamiltonian= b3lyp_G ] \
mm_theory=dl_poly : [ list \
cutoff= 1000 \
list_option=none \
exact_srf=yes \
use_pairlist=yes \
mxlist=30000 \
mxexcl=2000 \
scale14 = { 1.0 1.0 } \
atom_types= $types \
use_charmm_psf=yes \
charmm_psf_file= ${sys_name_id}.psf \
charmm_parameter_file= $prm \
charmm_mass_file= $top ] ]

#################################
#specification of atoms for scan#
#################################
# A: atom that hosts the displaced atom at the beginning
# B: displaced atom
set A 266
set B 4402

##################################
#scan parameters #
##################################
set stepnum 7

#set initdist [interatomic_distance coords= sam_0.c i=$A j=$B unit=angstrom ]
set initdist 2.4
set incr -0.2

for {set i 0} { $i < $stepnum} {incr i} {

set stepper [expr ($initdist + $incr*$i)*1.89036 ]

#coordinates from last minimization
set atomA [ get_atom_entry atom_number=$A coords= ${sys_name_id}_$i.c]
set atomB [ get_atom_entry atom_number=$B coords= ${sys_name_id}_$i.c]
#set atomA [ get_atom_entry atom_number=$A coords= scan_$i.c]
#set atomB [ get_atom_entry atom_number=$B coords= scan_$i.c]
set symbol [ lindex $atomB 0]

#create bond vector
set xa [ lindex $atomA 1]
set xb [ lindex $atomB 1]
set ya [ lindex $atomA 2]
set yb [ lindex $atomB 2]
set za [ lindex $atomA 3]
set zb [ lindex $atomB 3]

set xx [ expr $xb - $xa ]
set yy [ expr $yb - $ya ]
set zz [ expr $zb - $za ]

set length [ expr sqrt($xx*$xx + $yy*$yy + $zz*$zz) ]

set xxx [ expr $xx / $length ]
set yyy [ expr $yy / $length ]
set zzz [ expr $zz / $length ]

#displace atom along bond vector
set newx [ expr $xa + $xxx * $stepper ]
set newy [ expr $ya + $yyy * $stepper ]
set newz [ expr $za + $zzz * $stepper ]
set newat [subst {$symbol $newx $newy $newz}]
replace_atom_entry coords=${sys_name_id}_$i.c atom_number=$B atom_entry= $newat

#-------------------------------------------------------------
# optimize geometry with distance A-B fixed
#
dl-find coords= ${sys_name_id}_${i}.c \
result= ${sys_name_id}_[expr ($i+1)].c \
coordinates= cartesian \
optimiser= lbfgs \
tolerance= 0.0009 \
tolerance_e= 0.00001 \
residues= ${pdbresidues} \
residues= ${myresidues} \
constraints= {{bond 226 4402}} \
active_atoms= $active_atom_numbers \
list_option= full \
maxcycle= 999 \
maxene= 10000 \
maxstep= 0.5 \
dump= 100 \
theory= hybrid : "$hybrid_args"

if { $i == 0 } {
set noofqmatoms [get_number_of_atoms coords= hybrid.${qm_theory}.coords]
puts $control_input_settings "No. of atoms in the QM region: $noofqmatoms"
if { $incr < 0 } {
puts $control_input_settings "The scanned parameter is decreasing in value."
} else {
puts $control_input_settings "The scanned parameter is increasing in value."
}
}
flush $control_input_settings

puts $control_input_settings "======================================"
puts $control_input_settings "Structure [expr ($i+1)]"
# R1 = scanned distance
set R1 [interatomic_distance coords= ${sys_name_id}_[expr ($i+1)].c i=$A j=$B unit=angstrom ]
set R1 [expr $R1]
puts $control_input_settings "Scanned parameter (Distance $A-$B): $R1"

exec cp alpha alpha_$R1
exec cp beta beta_$R1
exec cp energy energy_$R1
exec cp gradient gradient_$R1
exec cp coord coord_$R1
exec cp natural natural_$R1

#Write pdb file
read_pdb file= ${sys_name_id}.pdb coords=dummy.coords
write_pdb file= ${sys_name_id}_$R1.pdb coords= ${sys_name_id}_[expr ($i+1)].c
write_xyz file= ${sys_name_id}_qm_$R1.xyz coords=hybrid.${qm_theory}.coords

set qmenergy [exec tail -2 energy | head -1 | cut -c7-22]
set qmenergy [expr $qmenergy]
puts $control_input_settings "Total QM Energy: $qmenergy"
set qmmmenergy [ get_matrix_element matrix=dl-find.energy indices= {0 0 } ]
set qmmmenergy [expr $qmmmenergy]
puts $control_input_settings "Total QM/MM Energy: $qmmmenergy"

flush $control_input_settings

#---------------------------------------------------------------------------

}

# cleanup
catch {file delete dummy.coords}
flush $control_input_settings
close $control_input_settings

exit

tomkeal
tomkeal's picture
Scan job error at cartesian coordinates

Hi,

If you are using Cartesian coordinates rather than HDLCs, you won't be able to apply a constraint on the bond distance with the "constraints=" option as these are specifically HDLC constraints. Instead, you can use the "restraints=" option in DL-FIND to apply a penalty term to the relevant bond distance. See here for how to define this:

https://www.chemshell.org/sites/www.chemshell.org/files/docs/tcl-chemshe...

Best wishes,

Tom

Log in or register to post comments