Scan job showing residue conversion error

4 posts / 0 new
Last post
Scan job showing residue conversion error

Does anything know anything about "residue conversion error" in chemshell ?
I am trying to do a scan between two atoms 147 & 495 which lie in residue number 5 and 16 respectively.
After 5-6 steps of scan job my log file displays this error:

Cyclic failure at residue 4, stopping

Residue conversion error
My INPUT file looks likes this:

# specify the name of the system

exec mkdir scan_job_0
set location ./scan_job_0

set control_input_settings [open ${location}/SUMMARY.txt a]
puts $control_input_settings "Summary of scan."

global sys_name_id
set dir .
set sys_name_id rna

exec cp rna.c ${location}/scan_0.c

source ${dir}/hybrid2.tcl
source ${dir}/orca.tcl
#load the connectivity
load_connect_from_psf ${dir}/${sys_name_id}.c ${dir}/rna_solvated.psf
#define force field
set top { top_all36_na.rtf toppar_water_ions.rtf }
set prm ./par_mod.inp
source ${dir}/save_${sys_name_id}.chm

source qmatoms
source act
source pdbresidues

set myresidues [ inlist function=combine residues= $pdbresidues sets= {RNAA5 RNAA16} target=fatone ]
#QM/MM energy
# Optimisation
global qm_theory
set qm_theory orca

fragment hybrid.${qm_theory}.coords new persistent
set hybrid_args [ list \
coupling= shift \
groups = $groups \
cutoff=1000.0 \
atom_charges= $charges \
qm_theory=orca : [ list \
hamiltonian=hf \
executable=/global/orca_4_0_1_2_linux_x86-64_openmpi202/orca \
nproc=12 \
maxcyc=20 ] \
qm_region= $qmatoms \
debug=no \
mm_theory=dl_poly : [ list \
list_option=none \
debug=no \
exact_srf=yes \
use_pairlist=no \
mxlist=16000 \
cutoff=1000 \
mxexcl=2000 \
debug_memory=no \
scale14 = { 1.0 1.0 } \
atom_types= $types \
use_charmm_psf=yes \
charmm_psf_file= ${dir}/rna_solvated.psf \
charmm_parameter_file= $prm \
charmm_mass_file= $top ] ] \

#optimize geometry

set A 147
set B 495

#scan parameters
set stepnum 20

set initdist [interatomic_distance coords= ${location}/scan_0.c i=147 j=495 unit=angstrom ]

set incr -0.1

for {set i 0} { $i < $stepnum} {incr i} {
#1.89036 changes angstrom to bohr
set stepper [expr ($initdist + $incr*$i)*1.89036 ]

set atomA [ get_atom_entry atom_number=$A coords= ${location}/scan_$i.c]
set atomB [ get_atom_entry atom_number=$B coords= ${location}/scan_$i.c]
set symbol [ lindex $atomB 0]

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=${location}/scan_$i.c atom_number=$B atom_entry= $newat

# optimize geometry with distance A-B fixed

dl-find coords= ${location}/scan_${i}.c \
result= ${location}/scan_[expr ($i+1)].c \
maxcycle=20 \
residues= $myresidues \
coordinates= hdlc \
active_atoms= $act \
maxene=10000 \
dump=100 \
list_option=none \
optimiser=lbfgs \
maxstep=0.5 \
tolerance=0.0009 \
tolerance_e=0.00001 \
constraints= {{bond 147 495}} \
theory= hybrid : "$hybrid_args"

#save output
exec cp ${location}/scan_[expr ($i+1)].c ${location}/scan_[expr ($i+1)]_0.c
# write summary to file
puts $control_input_settings "======================================"
puts $control_input_settings "structure [expr ($i+1)]"

set energy [ get_matrix_element matrix= indices= { 0 0 } ]
puts $control_input_settings [format "Energy:%14.6f" $energy]

set R1 [interatomic_distance coords= ${location}/scan_[expr ($i+1)].c i=$A j=$B unit=angstrom ]
puts $control_input_settings [format "Distance R1(A-B):%4.3f" $R1]

flush $control_input_settings


# cleanup
catch {delete_object hybrid.orca.coords}
catch {file delete dummy.coords}
catch {file delete hybrid.orca.coords}
catch {file delete conn.temp}
flush $control_input_settings
close $control_input_settings


It would be very helpful if someone can help me in this.
Thank you.

tomkeal's picture


Without knowing what is in residue 4 it's not obvious what is going on, so it would worth checking what that residue consists of and visualising its trajectory to see if anything strange is happening. If it's not obvious why the conversion is failing, you could alternatively optimise in Cartesians, applying a restraint to maintain the interatomic distance in place of the HDLC constraint.


Problem when using Cartesians

When I used Cartesian coordinates my interatomic distance (between A and B) does not get updated correctly , so that's why I used hdlc coordinates and the distance got updated correctly. Also I do not have a bond between between atom A and atom B , does the restraint option still work ? or should I use constraints only ?
Okay I will try visualising the trajectory for this .

Thank you so much for replying.

tomkeal's picture
Restraints are completely

Restraints are completely general as they are implemented as an energy penalty term, so they do not require a bond to be defined between the atoms in question, or the use of any particular coordinate system. They are less precise than HDLC constraints (which fix the distance at a given value), but should be good enough for the purpose of the scan.


Log in or register to post comments