ComQum

ComQum
Comqum is a method to do combined QM/MM calculations.
In principle, it is a general procedure for any QM and MM programs. In practice, it is currently implemented for Amber and Turbomole.

The proper references of ComQum are:

U. Ryde (1996) "The coordination of the catalytic zinc ion in alcohol dehydrogenase studied by combined quantum chemical and molecular mechanical calculations". J. Comput.-Aided Mol. Design 10, 153-164.
U. Ryde & M. H. M. Olsson (2001) "Structure, strain, and reorganization energy of blue copper models in the protein". Intern. J. Quant. Chem., 81, 335-347.

The current version is ComQum-01.
The interface programs are located in:

/home/bio/Bin/Gfortran

The code is in /home/bio/Comqum.

ComQum is a free interface program that you can get from Ryde by request.
However, note that to run ComQum, you also need the softwares Turbomole and Amber which are commercial.

The only available documentation is this file, the ComQum page and a detailed specification of the ComQum file formats: comqum01.pdf.

How to start a ComQum calculation (new Amber >6 version)
Form to fill in

Start with an equilibrated crystal structure (with hydrogen atoms) stored in pdb format with charges (follow the instructions on this page).

So far, we have always run ComQum with sphaerical systems without periodic boundary conditions.
In principle, it should probably be possible to run with periodic boundary conditions, but this is full of technical problems.
You must start with pdb structure with the QM system in the centre of the periodic box and all molecules wrapped into the central box (otherwise the point charges will be strange).
Check also that tleap does not move the atoms (run step 8 on the input pdb file).
After this, the set-up is normally without problems, but you need to change sander.in2 and 3 to employ periodic boundary conditions and reduce the cutoff to ~12 Å:
However, sander often gives infinite energies, owing to problems with the wrapping.
It is also somewhat unsatisfactorily that the MM calculations are periodic with a finite cutoff, but not the QM calculations.

If you already have a periodic octahedral simulation, you can truncate it to a spherical system:
cpptraj prmtop ptraj.in
with the input file (ptraj.in):
trajin mdrest3 trajout mdrest3-wrap restart image origin center familiar com :98 go

You can get the central residue (98 in this example) by running changepdb, command cap on the pdbout file. First remove all water molecules with command rw and a radius of 0 Å.

changeparm <<EOF
prmtop
p
pdb3
m
mdrest3-wrap
q
EOF

and then cut it to a spherical system with the cut command of changepdb (use cap first to decide a proper radius, you should include the whole protein and get a spherical system of water molecules, typically 10-15 Å smaller than the cap radius)
changepdb
pdb3
m
p
cu

0 0 0
45
n
pdb4-cut
n
Change the name of metal-bound water molecules to WAM and move them before the other water molecules. Normally already done
When selecting the quantum system, note that if there is a hydrogen bond from the back-bone NH or CO groups of a junction residue to the quantum system, this group should be included in the quantum system (because they otherwise get strange charges).
Normally not needed: Ensure that you have the desired junctions in the file juncfactor
(cp /home/bio/Data/Comqum/junctfactor .)

For each junction, it should contain the following entries:
Residue name, junction atom, junction neighbour, desired Amber atom type of junction atom;
Next row: jtype (indicating method and basis set) and ideal bond lenght with this method
The format is free, but fields must be delimited by spaces.

$version
2006

ASP CB   CA   HC # Asp OOC-CH3 junction, ff94
7   1.1065d0      # Acetate, BPRi/6-31G*, Cs, C-H in plane

Note that the entries depend on the Amber force field to use (the desired atom type) and the QM method to use (the ideal bond length).

Obsolete: use instead the file /home/bio/Data/Comqum/junctfactor_old
ASP CB   CA   ACA 1.5260d0 # CT-CT in parm99.dat
7   1.1065d0               # Acetate, BPRi/6-31G*, Cs, C-H in plane

Note that this file is specific for a certain force field (e.g. Amber ff99), although this is not explicit.
Run pdbtocomqum

Select a logfile
Enter the name of the pdb file
Do not search for short contatcts
Enter a new title (can be done in a file, see next point).
Select the quantum system, including the junction atoms that will be changed to H atoms (you first give the number of the amino acid, then the number of the atoms in this amino acid; one <CR> let you select a new amino acid; another <CR> ends the selection).

sample file

Select the radius for system 2 (MM system that is optimised). Typically ~6 Å.
Select the radius for system 3 (MM system that is fixed). Typically the rest of the protein.
Select amber (default)
Select the basis set (i.e. the type of junction, typically 12 (TPSS/def2-SV(P)).
Use the new format (default).
Enter the junction atoms (those that will become hydrogen atoms in the quantum system). This can also be done in the file mentioned in point 5.
Do not remove the charge of neighbours (answer n).
Do not redistribute any deleted charges (answer n).

If you nevertheless decided to remove charges of the neighbours:
Check that the net charge outside system 1 is correct (for normal side-chains, it should be 0, but for backbone fragments and strange groups, it may be <>0).
Check that the change in the charges is not too large. If so, remove more neighbours. For example, you may want to remove charges of hydrogens if you have already removed the charge of their parent heavy atom.

Run comqumtoturbo and comqumtoamber

comqum.dat comqum.qcout coord logfil pointcharges sander.in2 comqum.mmout control leap.in pdb.in3 sander.in1 sander.in3

Check the file coord - it often gives problems with metals (not zn, fe, ni).
Edit the file leap.in to insert unusual residues (i.e. insert extra addPrep and addAmberParm statements) and check the force field. Sample leap.in file.
Normally, it is best to use the same leap.in file as used for the equilibration run, but you must remove the solvateCap command, change the input pdb file to pdb.in3, and change the output files to prmtop3 and prmcrd3:
x=loadpdb pdb.in3saveamberparm x prmtop3 prmcrd3
Then, run tleap:

tleap -f leap.in

Check that you do not get any errors or warnings, especially after the commands
x=loadpdb pdb.in3saveamberparm x prmtop3 prmcrd3
Many errors come from missing or extra TER rows in the pdb file - leap will connect all consecutive residues unless they are separated by a TER row.

If you use non-standard charges for any residue, change the charges of system 1-3 by changeparm command m on prmtop3, using a pdb file with the correct charges (i.e. the file used as input to pdbtocomqum). Note that leap will change back to standard charges, even if the charges were modified in the input pdb file.

This is done by command cqcharge, i.e.
changeparm <<EOF
prmtop3
m
comqum.pdb
w

q
EOF
Normally not needed: Add the desired atom types of the junction atoms into file comqum.dat (if missing)
after the bond length.
Sample comqum.dat file.
Unusual: If there is no junctions in the QM system, run command
cqtrunc_no_junc
Run cqprep and follow the instructions. You will do the following:

Run cqtrunc to generate the file prmtop1 from file prmtop3 using the command
Remove charges from prmtop1 and write out pdbout1 with changeparm (done automatically)
Remove charges of system 1 from prmtop2 and write out pdbout3 with changeparm (done automatically).
Insert a cap into prmtop3 with changeparm - use the data from the equilibration. It is very important that you did not have a cap in prmtop3 before running cqprep.
Run a test of sander.in1 (write down the number of atoms and the energy). Continue the run by writing a q
Run a test of sander.in2 (this may take several minutes; write down the number of atoms and the energy). Continue the run by writing a q

Run checkcoord to see that the correspondance list works (check that there are no warnings (Amber may change the order of the atoms; if so, you have to reorder them also in the original pdb file and rerun everyting from pdbtocomqum (point 3). There should be no warnings. If you get warnings forr HB-atoms, e.g. in CYM, run fixcoord1 (fixcoord1_amberin; fixcoord1_turboin; fixcoord1; fixcoord1_amberout) and see if the program gives errors. If so, you have to swap coordinates in the prmcrd1 file.
Run a test of sander.in3 (write down the number of atoms and the energy

If you have two junction atoms that are were bound in the full (untruncated system), you need to delete the corresponding bond in the prmtop1 file. This is done with changeparm, command bb on file prmtop1.

Normally not needed: If you run several states and want to compare the energies it is absolutely mandatory that you have exactly the same system2 (i.e. the same belly in sander.in2) in all calculations. You do that in pdbtocomqum by adding or removing residues to system 2 (step 6).
However, it often enough to copy the same sander.in2 and sander.in3 files to all directories (check with diff that they are identical).
Run define to set up the QM calculation (more instructions; explicit commands):

Make an UFF hessian - coordinate menue, command ff, set the charge and then hit return; when exiting the coordinate menu, you need to explicitly say n, you do not want to define internal coordinates;
Define the basis set.
Set up the wavefunction (eht)
Optionally change the DFT functional
Turn on RI
Turn on dispersion with bj damping.

Remove unneccesary files (an make a backup) by command
cqzip
(gzip *; cqgunzip; mkdir Gz; mv *.gz Gz; cp * Gz; gzip Gz/*&)

The following 15 files are left (for UHF jobs, there are 16 files and alpha and beta replace mos; perhaps also auxbasis):
basis control mos prmcrd1 prmtop1 prmtop3 sander.in2 sander.out3
comqum.dat coord pointcharges prmcrd3 prmtop2 sander.in1 sander.in3
Typically you first run a set of protein-fixed calculations (default input; i.e.
$protein fixed
is in the file comqum.dat file.
Then run
comqum -c 200
or (with the RI approximation):
comqum -ri -c 200Always check the convergence of your results. We use the standard Turbomole quite floppy convergence criteria, which means that you sometimes can get random convergence. Therefore, check that the energy really is converged.The final results are in the prmcrd3 and in the Turbomole files.
You can build pdb files from those (pdb and pdb3) using the command
cqtopdb
Then, you may run protein-free (i.e. you relax system 2 with MM). Set (in a new directory)
$protein free
in file comqum.dat.

If you have a residue with net charged groups outside the QM system, you need to insert into comqum.dat (change the residue number and the desired charge).
$residue_charge_corr
1
652 -1

Then, you can start the program with
comqum -ri -backup -c 200

Sometimes, the calculation crashes, saying that the charge of the system has changed.
This normally means that your leap charges (typically of the QM system) were not correct.
Check what is the correct charge of the protein and if it is not correct in the prmtop3 file, use cqtopdb to write a pdb file with the prmtop charges, then change the charge of one QM atom to get the correct total charge of the protein, and insert those charges into the prmtop3 file with changeparm, command m.

You can test this part of comqum by running:
ridft >logd
sed -i 's/$point_charges/#point_charges/' Str/control
ridft -proper >mulliken
sed -i 's/#point_charges/$point_charges/' Str/control
fixcharge_turboin>>fixc
fixcharge_amberin>>fixc
fixcharge>>fixc
fixcharge_amberout>>fixc

Sometimes it can be quite a pain to find errors in the charges. Here are some hints.
- The charge of all atoms in residues outside the QM system must be an integer.
- Charges of residues with junctions must be determined for the same size of the QM system (or smaller; dummy charges are OK in the QM system, but net charge must be right).
- In the output of fixcharge:
- Sum of change must be 0 (=oldch + diff - QMch)
- Sum of QMch must be integer
- The other three sums can be non integer, but newch=oldch and diff+newch=QMch
- For each residue:
- diff=(RestSum-Corr)/#junc
- Sum of oldch + RestSum must be an integer
- OldSum and Corr shall sum up to an integer (not always)
For big proteins, it is not possible to run with an infinite cutoff (cut=1000.0 in sander.in3). Then, set it to the largest possible number. It is then normally also necessary to run with more floppy convergence thresholds:
comqum -ri -energy 4 -gcart 2 -c 200and after convergence swith back to protein fixed and rerun comqum.
Always check the convergence of your results. We use the standard Turbomole quite floppy convergence criteria, which means that you sometimes can get random convergence. Therefore, check that the energy really is converged. This is especially important with protein free.
You should typically run qmmmpbsa on all ComQum results (it takes only ~1 h).
This is done by
qmmmpbsa-vac >logv
qmmmpbsa-ptch>logp
Note that you need the files, so you cannot copy the calculation to the temp disk of the batch computer.

This is a first and minimal step to make the QM/MM eneriges more stable (and to include a proper solvation).
Typically, you should also run calculations forth and back between the start and end states, especially if you have protein free, to ensure that the energies are stable (that the two states reside in the same local minima) or start QM/MM optimisations on several snapshots from a MD simulation.

Alternatively, you should consider to run QTCP or other QM/MM free-energy perturbation methods.

If you get the problem
Too large difference - use $residue_charge_corr
it often means that you have a different (incorrect) charge of the QM system in the prmtop3 file (for a protein-fixed calculation, this does not matter and is not noticed by comqum until you run qmmmpbsa).
Add the missing charge to one atom in the QM system in the pdb3 file (by hand) and then run:

changeparm <<EOF prmtop3 mpdb3 w prmtop3 q EOF

Then, you can run qmmmpdbsa as usual.
You can add restraints on a distance between two atoms using in comqum.dat (note the unusual units):
$restraints
atom1(1) atom2(1) desired_distance_in_A   force_constant_in_au
$restraints
1 3 1.00 1.0
Note that the atoms should be numbered according to the quantum system and not according to the total (MM) system.
A typical force constant is 1.0.
You may also use a higher force constant, e.g. 10.0 a.u. (which gives a smaller deviation from the target, typically 0.05 pm), but the you often get oscillations in the energies or slow convergence.

You can also restrain the difference between two distances (r_AB - r_BC):
$restraints
   atom1 atom2 atom3 desired_distance_in_A   force_constant_in_au
$restraints
1 2 3 1.00 1.0

You can also restrain an angle:
$restraints
   atom1 atom2 -atom3 desired_angle_in_deg force_constant_in_au
Note the minus sign before atom3 (to discern it from the previous restraint).
$restraints
1 2 -3 1.00 1.0
A force constant of 1-10 au seems to be OK.

You can also restrain a torsion angle:
$restraints
   atom1 atom2 atom3 atom4 desired_angle_in_deg force_constant_in_au
$restraints
1 2 3 4 1.00 1.0
Again, a force constant of 1 seems to be OK.
Similarly offset forces are added by (also in comqum.dat):
$offset forces
   atom1(1) atom2(1) force(1)
   ...
   atom1(n) atom2(n) force(n)
The atoms should be given by numbers, forces in atomic units (the negative of the output in relax, possibly calculated with
$interconversion on
   gradient cartesian --> internal
this works also with redundant internal coordinates if you first define the desired coordinates)
Note that the atoms should be numbered according to the quantum system and not according to the total system.
In order to calculate strain energies (the first two lines can be done with command str):
mkdir Str
cp control coord basis energy mos alpha beta auxbasis Str
cd Strdscf >logd (or ridft >logd)
To evaluate the effect of the protein, you should use the following four data points: The QM/MM energy, E(QM/MM), the QM energy with point charges, E(QM+ptch), found in the fourth column in the energy file, the QM vacuum energy at the QM/MM geometry (optained by the str script or by qmmpbsa-vac, E(QM1), and the QM energy optimised in vacuum, E(QM). Then:
E(QM1)-E(QM) is the strain energy of the QM system, i.e. the effect of the change in geometry in the protein.
E(QM+ptch)-E(QM1) is the effect of the point charges in the protein
E(QM/MM)-E(QM+ptch) is the MM energy. With protein fixed, it is usually small. If not, you should run the reaction forth and back until the energy difference converges. It is typically an artifact of different local minima.

Typically, the ptch energy dominates, and it can be further understood by using changepdb, command im2, on the two pdb3 files (for two states in an reaction) together with the corresponding charges in Ptch/logd. It will perform a per-residue decomposition of this energy (using a MM approximation, so the difference is not exactly equal to E(QM+ptch)-E(QM1), but normally quite close.

The MM (=MM12-MM1) energy difference between two states can be understood by running
changeparm, on parmtop2 command qmd (QM/MM difference)
Read in the corresponding prmtop1 file (they should apply to both states)
prmcrd3 and prmcrd1 of the first state
prmcrd3 and prmcrd1 of the second state
Further improvement of the energy differences can be obtained with big-QM and QTCP calculations.

There is a potential currently unsolved problem in all ComQum calculations:
If junction atoms are involved in improper dihedrals, there will be a mismatch between MM1 and MM3, because the movement of the junction atom will make the improper dihedral angle not equal in MM1 and MM3. It is normally not an important problem. We have seen problems from it only for very distorted structures (caused by other problems).
It is solved by using additive QM/MM.

The ComQum energy is:
    E_QM/MM = E_QM1+ptch23 + E_MM123_no_q1 - E_MM1_no_q
In the energy file, you find the E_QM/MM and E_QM1+ptch23 in the second and fourth columns:
   628 -16161.8463086333 15231.6278196500 -15322.2901877800     -3.4144500185           E_QM/MM            SCF_kin      E_QM1+ptch23              E_MM_opt
The two MM energies are in the files calcforce.out1 (E_MM1_no_q) and calcforce.out2 (E_MM123_no_q1), which can be accessed by
grep 'Total energy               =' calcforce.out?

Explicit Turbomole define commands

define
yffc -2nobb all def2-SV(P)bl*eht-2ydftfunc tpss*rionm 2500*dspbj**

How to set up a system with all groups within a specified radius within an original QM system

This can be done with the PMISP tools:

Start from a ComQum calculation of the original QM system (comqum.dat and pdb3 files)
If you have CYM or other strange residues, change them to standard residue names.
Run changepdb on pdb3
Enter command sp
Enter i
Enter the original QM system as ranges (first and last atom; same if no range), with the second atom with a negative number if it is not the last one:
For example:
$correspondance_list
    33
   531   534   535   536   537   538   539   540   541 1229 1231 1232
1233 1234 1266 1269 1270 1271 1272 1273 1274 1275 1276 1324
1327 1328 1329 1330 1331 1332 1333 1334 1443
becomes:
531
-531
534
-541
1229
-1229
1231
-1234
1266
-1266
1269
-1276
1324
-1324
1327
-1334
1443
1443
Enter the desired radius
Enter return once
Enter auto
Enter restop
Enter return
x
pdb3.xyz
Finish with q
Finally enter
pmisp 1 pdb3.xyz split.inp - -
This generates a large number of files. The only you need is
systa
\rm f*xyz *input prmcrd* comqum_a* c*xyz comqum_b.dat comqum_b1.dat gtot.xyz tot.xyz systb g1.xyz

To avoid so much output, you can also make a file dummy.inp with "0 0" in the first line and an empty second line and then run
pmisp 1 pdb3.xyz split.inp - dummy.inp

If you get when issuing the auto command
Sorry, residue 62 is disconnected.
This cannot be handled by a standard restop file.
Please do something about this.
Quitting...
It typically means that a residues contains two fragments (backbone and sidechain). Try to connect them by adding heavy atoms (in point 8 above).

This is the command cqtrunc

changeparm <<EOF
prmtop3
tr
comqum.dat
n
y
y
prmtop1
y
m
prmcrd3
prmcrd1
q
EOF

In command cqtrunc_no_junct one of the two consecutive y rows is deleted

This is`cqprep`

# Run cqtrunc echo "Running cqtrunc" cqtrunc # Zero charges in prmtop1 # This is now done by command cqzero1 echo "Zero the charges in prmtop1" changeparm <<EOF prmtop1 p pdbout1 m prmcrd1 z w prmtop1 q EOF # Zero the charges of the first QM system in prmtop2 # This is now done by command cqzero3 echo "Zero the charges of syst 1 in prmtop2" changeparm <<EOF prmtop3 p pdbout3 m prmcrd3 1 w prmtop2 q EOF# Insert the cap into prmtop3 echo " " echo "Now you should add a cap to prmtop3 if appropriate" echo "Enter: prmtop3; c; y; enter the appropriate values; w; <CR>; q" echo " " changeparm # Test calcforce calcforce <<EOF >calcforce.out1 prmtop1 m prmcrd1 1000 mden1 force1 EOF lesscalcforce.out1calcforce <<EOF >calcforce.out2 prmtop2 m prmcrd3 1000 mden2 force2 EOF
less calcforce.out2

# Use define to include the proper basis sets and get a starting wave function.

# DFT functional is already defined echo "Define basis set and mos by define" define # Check the corresponance lists so that Amber has not renumbered the atoms checkcoord # Test sander 3 (must be killed) echo "sander.out3: This job must be killed" sander -O -i sander.in3 -o sander.out3 -r mdrest3 -e mden3 -c prmcrd3 -p prmtop3 more sander.out3

Sample comqum.dat file

$title
Title
$protein fixed
$junction= 6
$junctions

     1      2     1.09800000000000 H1 
$correspondance_list
    15
  1976  1978  1979  1980  1981  1982  1983  4865  4866  4867  4868  4869
  4870  4871  4872
$junction_atoms
     1     2    1.38979963570127
$end

Note that $junctions is added by pdbtocomqum, whereas $junction_atoms is added by changeparm, command tr.

This is a schematic overview of comqum

# Start with a scf step then loop through the following four steps

# Gradient step (gradient calculation)
grad >logg
cp gradient grad1
offsetf >fixf
fixforce_turboin>>fixf
fixforce_amberin>>fixf
fixforce>>fixf
fixforce_turboout>>fixf
mv gradient g1
mv grad1 gradient

# Relaxation step
relax >logr
fixcoord1_turboin >fix1
fixcoord1_amberin>>fix1
fixcoord1>>fix1
fixcoord1_amberout>>fix1

# Molmech step
# if protfree then
   sander -O -i sander.in3 -o sander.out3 -r mdrest3 -c prmcrd3 -p prmtop3    cp mdrest3 prmcrd3
   fixcoord2_turboin >fix2
   fixcoord2_amberin>>fix2
   fixcoord2>>fix2
   fixcoord2_turboout>>fix2
#endif

# Scf step (energy calculation)
calcforce <<EOF >calcforce.out1 prmtop1 m prmcrd1 1000 mden1 force1 EOF calcforce <<EOF >calcforce.out2 prmtop2 m prmcrd3 1000 mden2 force2 EOF dscf>logd
grep 'Self energy' logd > selfenergy
fixenergy_turboin >fixe
fixenergy_amberin>>fixe
fixenergy>>fixe
fixenergy_turboout>>fixe
#if protfree then
   moloch>mulliken
   fixcharge_turboin >fixc
   fixcharge_amberin>>fixc
   fixcharge>>fixc
   fixcharge_amberout>>fixc
#endif
checkconv

These are the programs you need (apart from Turbomole and Amber):
forcetograd
fixcoord1
fixcoord2
fixenergy
fixcharge
comqum

In addition, you may need some of our accessory programs:
pdbtocomqum
cqprep
espprep
mkchargecorr
checkcoord
cqgunzip
changeparm
changepdb

This is no longer needed (Sept. 2016)
In addition, you need a modified version of sander (sanderf) which writes out the forces.
Three files in Amber/src/sander have been modified (Amber 8.0):
cp md.h md.h_orig
cp force.f force.f_orig
cp mdread.f mdread.f_orig

md.h
8c8
< parameter (BC_MDI=55)
---
> parameter (BC_MDI=54)
28c28
<       idecomp,icnstph,ntcnstph,maxdup,iwforc     !55
---
>       idecomp,icnstph,ntcnstph,maxdup     !54
42c42
<       ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup,iwforc
---
>       ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup

force.f:
five lines added after line 752 (1073 in Amber 9; 550 in Amber 7; 677 in Amber 5; almost at the end);
after the row:
call trace_exit( 'force' )

!------- This part added for ComQum --- ! Write out the forces if (iwforc.eq.1) then write(77,'(3e22.14)')(f(i),i=1,3*natom) write(77,*) endif !----------------------------------- return
end subroutine force

mdread.f
Three lines as this diff comparison shows:

84c84 < mmtsb_switch, mmtsb_iterations,rdt,icnstph,solvph,ntcnstph,iwforc & --- > mmtsb_switch, mmtsb_iterations,rdt,icnstph,solvph,ntcnstph & 266,268d265 (Before line ! Check to see if "cntrl" namelist has been defined) < ! Write forces? < iwforc=0 <

This is an example of a file with the quantum system (syst1)
New title
531
534-541
1229
1231-1234
1266
1269-1276
1324
1327-1334
1443531
1229
1266
1324

The format is :
A new title line
Atoms of the QM system (one atom or a range on each line)
Blank line
Junction atoms (one on each line)

This is an example of a Gaussian file for calculation of the fitted charges
%Chk=impc42cm.chk
#P HF/Sto-3G SCF(Vshift=1000) Charge=Bohr

Pc.ox/Cu2, ComQum on ImPc42CM, 16/24 A, DZpdf/6-31G*, 3/12-99
pcox.new, Plastocyanin, new equilibration with AMBER, 24 A CAP, 960624

    1    2
h      -2.34543400328041    7.56091572587130   -4.16389870653224
c      -2.51199881855881    8.35999606813360   -3.48399836140879
n      -1.51299928840743    9.17099568670494   -2.99999858904316
c      -2.10299901091925   10.03399528081969   -2.17699897611565
h      -1.58799925313351   10.80799491679282   -1.62799923432075
n      -3.43899838257314    9.80999538617113   -2.10099901185989
h      -4.09199807545487   10.32199514536783   -1.52499928276361
c      -3.72299824900256    8.74799588564985   -2.92699862337644
h      -4.69299779279318    8.29799609729338   -3.07899855188796
h       3.04952150153140    5.77052939610116   -3.07577470019319
c       2.27499893002440    6.42199697961172   -2.68199873860458
h       1.50799929075903    5.80699726885787   -2.21199895965449
h       2.71799872167310    7.06199667860760   -1.91899909745794
s       1.53599927759010    7.43899650129735   -3.92699815305749
h       2.91307908205073   10.89635139119452   -2.44756015205704
c       2.16999897940789   11.41399463177954   -3.03699857164136
n       0.90199957577231   10.89799487446412   -3.20299849356841
c       0.26899987348420   11.77099446387568   -3.94899814271048
h      -0.74299965055302   11.64499452313586   -4.30499797527693
n       1.03399951369021   12.83399396392664   -4.23299800913990
h       0.77999963315122   13.58899360883583   -4.85399771707183
c       2.25999893707918   12.64899405093564   -3.64899828380616
h       3.10999853730807   13.31499373770322   -3.68499826687468
h       2.51873747856117    8.39759526408039    0.78586489906895
c       1.59999924748968    7.82499631975424    0.73999965196398
h       1.82399914213824    6.83299678631063    0.34799983632901
h       1.20799943185471    7.70999637384092    1.75099917647152
s       0.34399983821028    8.58399596278216   -0.25599987959835
c      -0.87399958894124    7.21399660711912   -0.22199989558919
h      -0.42999979776285    6.29999703699063   -0.61599971028353
h      -1.18699944173141    7.02499669600940    0.80499962139325
h      -1.74799917788248    7.47799648295492   -0.81699961574942
cu      0.37499982363039    8.91299580804723   -3.05899856129434