Resp fitting [Cu(isaepy)H2O]+

  babel -ig09 Cuisaepy_RESP.log -omol2 Cuisaepy_RESP.mol2
1 molecule converted
1 info messages 25 audit log messages 
  e Cuisaepy_RESP.mol2 
  atomtype -i Cuisaepy_RESP.mol2 -o Cuisaepy.ac -f mol2 -p gaff

WARNING: atom type of     C (cc) and     N (nc) may be wrong
  espgen -i Cuisaepy_RESP.log -o  Cuisaepy.esp
  respgen -i  Cuisaepy.ac -o  Cuisaepy.respin1 -f resp1
  respgen -i Cuisaepy.ac -o Cuisaepy.respin2 -f resp2
  resp -O -i Cuisaepy.respin1 -o Cuisaepy.respout1 -e Cuisaepy.esp -t qout_stage1
  resp -O -i Cuisaepy.respin2 -o Cuisaepy.respout2 -e Cuisaepy.esp -q qout_stage1 -t qout_stage2
  antechamber -i Cuisaepy.ac -fi ac -o Cuisaepy.ac -fo ac -c rc -cf qout_stage2

For atom[3]:Cu, the best APS is not zero, bonds involved by this atom are frozen
Info: Bond types are assigned for valence state 3 with penalty of 2

WARNING: atom type of   C13 (cc) and    N2 (nc) may be wrong

  prepgen -i Cuisaepy.ac -o Cuisaepy_resp.prep -f int -rn LIG -rf LIG.res
  parmchk -i Cuisaepy_resp.prep -f prepi -o Cuisaepy.frcmod
Segmentation fault
  cat Cuisaepy_RESP.mol2 | awk '{ sum+=$9} END {print sum}'
0.9999

# Este comando grava as cargas ESP no arquivo mol2
babel -ig09 Cuisaepy_RESP.log -omol2 Cuisaepy_RESP.mol2

atomtype -i Cuisaepy_RESP.mol2 -o Cuisaepy.ac -f mol2 -p gaff

espgen -i Cuisaepy_RESP.log -o  Cuisaepy.esp

respgen -i  Cuisaepy.ac -o  Cuisaepy.respin1 -f resp1
respgen -i Cuisaepy.ac -o Cuisaepy.respin2 -f resp2
resp -O -i Cuisaepy.respin1 -o Cuisaepy.respout1 -e Cuisaepy.esp -t qout_stage1
resp -O -i Cuisaepy.respin2 -o Cuisaepy.respout2 -e Cuisaepy.esp -q qout_stage1 -t qout_stage2


antechamber -i Cuisaepy.ac -fi ac -o Cuisaepy_resp.ac -fo ac -c rc -cf qout_stage2
prepgen -i Cuisaepy_resp.ac -o Cuisaepy_resp.prep -f int -rn LIG -rf LIG.res


parmchk -i Cuisaepy_resp.prep -f prepi -o Cuisaepy.frcmod

cat Cuisaepy_RESP.mol2 | awk '{ sum+=$9} END {print sum}'

Resp fitting, leap - Cuisaepy

Resp Fitting:

Com o resultado das cargas resp do gaussian deve-se fazer:

# Converter a estrutura g09 em mol2 (observe qual carga foi adicionada no mol2: Mulliken ou ESP)
babel -i g09 cuisa_opt.log -o mol2 cuisa_opt1.mol2

# Conferir as cargas
cat  cuisa_opt1.mol2 | awk '{ sum+=$9} END {print sum}'

#Modificar o mol2 para ser reconhecido pelo gaff (Cu)
sed -i 's/CU /Cu /g' cuisa_opt1.mol2

#Generate ac file with atomtypes (not correctly yet);
atomtype -i cuisa_opt1.mol2 -o cuisa_opt.ac -f mol2 -p gaff

## Mensagem no output:
## WARNING: atom type of     C (cc) and     N (nc) may be wrong


#Take ESP charges from the g09 output file
espgen -i cuisa_opt.log -o cuisa_opt.esp

#Fitting procedure
respgen -i cuisa_opt.ac -o cuisa_opt.respin1 -f resp1

respgen -i cuisa_opt.ac -o cuisa_opt.respin2 -f resp2

resp -O -i cuisa_opt.respin1 -o cuisa_opt.respout1 -e cuisa_opt.esp -t qout_stage1

resp -O -i cuisa_opt.respin2 -o cuisa_opt.respout2 -e cuisa_opt.esp -q qout_stage1 -t qout_stage2

#Use antechamber to make prep file
antechamber -i cuisa_opt.ac -fi ac -o cuisa_opt_resp.ac -fo ac -c rc -cf qout_stage2

#ERRO (ver)

prepgen -i cuisa_opt_resp.ac -o cuisa_opt_resp.prep -f int -rn CIE -rf CIE.res

## Eu apaguei as ultimas 3 linhas com X  do cuisa_opt_resp.prep (acho que é o H2O....)
#cp cuisa_opt_resp.prep cuisa_opt_resp_mod.prep


# Obter o arquivo frcmod e neste ponto o frcmod não deve reconhecer o Cu!!!!!!!
parmchk -i cuisa_opt_resp.prep -f prepi -o cie.frcmod



cat  cuisa_opt_resp.prep | awk '{ sum+=$11} END {print sum}'

Aqui deve-se colocar os parâmetros obtidos pelo Force Matching

No arquivo frcmod eu coloquei:
MASS
Cu 63.6

NONB
  Cu          1.2000  0.0500      (está certo????)

AMANHA:
/local/GRS/Users/brown/complexes/Cuisaepy/CuisaepyH2O+/1/RESP/New

1) Ver no xleap quem sao os atom types e colocar os valores das distancias e angulos e diedrais
2) Acertar a carga que era a do H2O (pode distribuir entre o Cu e os sitios ligantes......
3) Rodar uma minimizacao e ver se o complexo é estavel
4) Ver se é possivel que haja a troca entre cotovelos.....



5) /local/GRS/Users/brown/complexes/Cuisaepy/CuisaepyH2O+/1/RESP/IOP
Obter o resp diretamente pelo antechamber (ver pag 148 do manual)

###################
# topology files
###################

#copiar os arquivos cuisa_opt_resp.prep e cie.frcmod para a pasta 3-LEAP

#Run tleap and set some atom types
tleap -f leaprc.gaff
source leaprc.ff99SB

loadamberprep cuisa_opt_resp.prep

savemol2 CIE cie.mol2 1

saveoff CIE cie.lib
loadamberparams cie.frcmod
saveamberparm CIE cie.top cie.crd
savepdb CIE cie.pdb
addions CIE Cl- 0
solvatebox CIE TIP3PBOX 14
saveamberparm CIE cie_solv.top cie_solv.crd
savepdb CIE cie_solv.pdb
quit

mpirun -np 8 pmemd.MPI -O -i min1.in -o min1.out -p cie_solv.top -c cie_solv.crd -r cie_solv_min1.rst -ref cie_solv.crd &

mpirun -np 8 pmemd.MPI -O -i min2.in -o min2.out -p cie_solv.top -c cie_solv_min1.rst -r cie_solv_min2.rst &

mpirun -np 8 pmemd.MPI -O -i md1.in -o md1.out -p cie_solv.top -c cie_solv_min2.rst -r cie_solv_md1.rst -x cie_solv_md1.mdcrd -ref cie_solv_min2.rst &

mpirun -np 8 pmemd.MPI -O -i md2.in -o md2.out -p cie_solv.top -c cie_solv_md1.rst -r cie_solv_md2.rst -x cie_solv_md2.mdcrd &

Resultados Parametrização [Cu(isaepy)]+

First coordination shell of [Cu(isaepy)]+
water and name O and same residue as (within 3.15 of name Cu)

Energia ptencial da simulação do Cuisaepy+

g(r) do cobre e do oxigênio das moléculas de H2O. O primeiro pico está em 2.15A com ppulação de 2.75. O primeiro vale está a 3.15A. A segunda camada de solvatação está a 4.55A com população de 0.9.

Cuenim - simulation results

Dinâmica molecular clássica do [Cu(enim)H2O]2+

Primeira e segunda camadas de solvatação através de dinâmica clássica 

Dinâmica molecular quântica do [Cu(enim)H2O]2+

resname SOLV and same residue as (within 2.15 of resname LIG)

Distâncias médias do cobre aos oxigênios da água (QMMM)

Cu_solv_inferior_dist.dat 
Media (Desvio): 2.14  0.08
Cu_solv_superior_dist.dat 
Media (Desvio): 2.16  0.10

Distância da H2O equatorial (QM): ~ 2.075 A

Distâncias Cu-OH2

[Cu2+(Cl-)2](H2O)2    2.15 (0.14)

[Cu(enim)H2O]2+        2.075 (equatorial - QM) 
axiais
2.16  0.10 (superior) 
2.14  0.08 (inferior)

[Cu(isaenim)]2+          
2.18 (superior) 
2.11 (inferior)

Copper complex parametrization

Empirical Force Fields for Biologically Active Divalent Metal Cations in Water
C. Satheesan Babu, and Carmay Lim
J. Phys. Chem. A 2006, 110, 691-699

Informações:
Raio de van der Waals para o Cu2+
Este método não descreve a distorção Jahn Teller,  o efeito do elétron desemparelhado do cobre alargando as distâncias axiais ou equatoriais.

vdW parameters for Cu2+

             e        Rmin/2      deldelG         CN     Cu-O
         (kcal/mol)  (angstron)  (kcal/mol)       #   (angstron)
Model A 0.0734   0.8685   -13.3 (-0.1)  6    1.90
Model B 0.1305   1.0763   -13.5 (-0.3)  6    2.03
Model C 0.0427   1.0330   -59.9 (+1.0)  6    1.94

Qual usar? Seria melhor utilizar o modelo que descreve melhor os parâmetros estruturais:
A escolha foi a do modelo C.

Preciso fazer a conversão para o vdW do Cu2+?? Já que este é para o sistema Cu-H2O?????

Observe que no arquivo parmtop uso o Rmin/2 como está listado na tabela acima


Discussion about CuH2O structural geometry
from J. Phys. Chem. A, Vol. 114, No. 35, 2010
Development and Validation of a ReaxFF Reactive Force Field for Cu Cation/Water

Interactions and Copper Metal/Metal Oxide/Metal Hydroxide Condensed Phases

For a long time, it was generally accepted that the copper(II) ion is 6-fold coordinated in aqueous solution with a tetragonally distorted octahedral structure, as would be predicted from the Jahn-Teller effect for a d9 electronic configuration. 

This was challenged in 2001 by Pasquarello et al. who proposed a 5-fold coordination based on first-principles MD and neutron diffraction.[10] Their simulation indicated frequent exchanges of square pyramidal and trigonal bipyramidal clusters with five equal Cu-O distances. One year later Persson et al. performed extensive EXAFS and large angle X-ray
scattering (LAXS) experiments and again suggested that the ion was 6-fold coordinated.[14] They modeled the copper-water cluster using several different geometrical configurations and showed that a tetragonally distorted octahedron gave the best fit to the experimental data. In 2003, a quantum-mechanical/molecularmechanical (QM/MM) study by Schwenk et al. also indicated a 6-fold coordination with Jahn-Teller distortion.[12,13] However, another Car-Parrinello molecular dynamics (CPMD) simulation performed by Amira et al.[4] and analysis of the XANES portion of the XAS spectra by Benfatto et al.[15] and Frank et al.[16] suggested that the hydrated copper(II) ion can be better represented as a five-coordinate square pyramidal structure with one elongated axial water molecule.

More recently, Chaboy et al. performed XANES and EXAFS spectroscopy of copper ion hydration [17,18] and argued that “neither the classical Jahn-Teller geometry nor other six- or fivefold coordinated geometries may be proposed unambiguously as the single preferred structure in solution”.[18] According to them, a dynamic view of different geometries exchanging with each other is the right description of copper ion hydration. A review of the structural parameters of Cu2+ aqueous complexes is referenced.[19]

 MASS 
> CU 63.55 
> 
> BOND 
> nb-CU 75.582 2.027 
> o-CU 84.879 1.907 
> 
> ANGLE 
> ha-ca-nb 52.970 116.000 same as ha-ca-n2 
> ca-nb-CU 50.523 114.913 from g09 calculations 
> 
> DIHE 
> ca-nb-cp-ca 1 4.800 180.000 2.000 
> ca-nb-cp-cp 1 4.800 180.000 2.000 
> nb-CU-nb-cp 1 14.800 180.000 2.000 
> 
> IMPROPER 
> o- o-CU-nb 1.1 180.0 2.0 
> 
> NONBON 
> CU 1.0330 0.0427 LJ parms from 
> http://dx.doi.org/10.1021/jp054177x 

Enim - amber

# Check the parameters not included

parmchk -i case.prep -f prepi -o case.frcmod
parmchk -i case.mol2 -f mol2 -o case.frcmod

# Check in mol2 file whether the total charge is the expected one:
 
cat case_RESP.mol2 | awk '{ sum+=$9} END {print sum}'
1.9996
# Add the remaining charge to be equal to 2.00


# Observe that gaff does not recognize Copper atom

more $AMBERHOME/dat/leap/parm/gaff.dat

# Check the parameters not included in the frcmod file

This command provides a frcmod file:

remark goes here
MASS
Cu 63.6                       ATTN, need revision


NONBON
  Cu          1.2000  0.0500  ATTN, need revision



Run leap


tleap -f leaprc.ff99SB
source leaprc.gaff

LIG = loadmol2 case.mol2
loadamberparams case.frcmod
saveoff LIG case.lib
saveamberparm LIG case.top case.inpcrd

addions LIG Cl- 0
saveoff LIG case_CI.lib
saveamberparm LIG case_CI.top case_CI.inpcrd

solvatebox LIG TIP3PBOX 14
saveoff LIG case_solv.lib
saveamberparm LIG case_solv.top case_solv.inpcrd



tleap -f leaprc.gaff
source leaprc.ff99SB
loadamberprep case.prep




5) Run sander

mpirun -np 8 sander.MPI -O -i min1.in -o min1.out -p case_solv.top -c case_solv.inpcrd -r case_solv_min1.rst -ref case_solv.inpcrd &


mpirun -np 8 sander.MPI -O -i min2.in -o min2.out -p case_solv.top -c case_solv_min1.rst -r case_solv_min2.rst &

mpirun -np 8 sander.MPI -O -i md1.in -o md1.out -p case_solv.top -c case_solv_min2.rst -r case_solv_md1.rst -x case_solv_md1.mdcrd -ref case_solv_min2.rst &


mpirun -np 8 sander.MPI -O -i md2.in -o md2.out -p case_solv.top -c case_solv_md1.rst -r case_solv_md2.rst -x case_solv_md2.mdcrd &

http://ambermd.org/antechamber/antechamber.html

DNA major groove intercalator

Major groove intercalator

Minor groove binder

Simulação Cu(H2O)

Teste realizado com dinâmica molecular clássica com os seguintes parâmetros não ligados de XXns
(500 frames gravados)

NONBON
  Cu          1.2000  0.05000             ATTN, need revision

Carga 2+

g(r) Cobre- Cloreto

g(r) Cobre - Oagua

Como a carga inserida para o cobre foi +2e, e os cloretos (-1), os cloretos permaneceram próximos ao cobre durante a simulação.

__________________________________________

[Cu2+(Cl-)2](H2O)x  - Parte quântica - Parte clássica



Teste realizado com o CPMD por xx ps com (5000 frames gravados).

TIMESTEP
 5.0
MAXSTEP

 10000
time step 1 a.u. = 0.02418884326505 fs
~ 1.2 ps

Cu_Cl_2.dat    Media (Desvio): 2.11 0.08

Cu_Cl_3.dat    Media (Desvio): 2.10 0.08

Cu_solv189.dat Media (Desvio): 2.16 0.14

Cu_solv297.dat Media (Desvio): 2.15 0.14

g(r) Copper+Ow

g(r) Copper+Cl-

Simulação DNA+Cuisaepy+

# Copiar  modelo que será utilizado no docking

tleap -f leaprc.gaff
source leaprc.ff99SB


# CARREGAR O COMPLEXO DE COBRE (o loadoff deve vir antes do loadpdb !!!!!!)
loadoff cie.lib  
CIE = loadpdb cie.pdb
loadamberparams cie.frcmod
check CIE


# CARREGAR O DNA
DNA = loadmol2 3GC.mol2
savepdb DNA DNA.pdb
check DNA

# COMBINAR OS DOIS RESIDUOS
CIED = combine { CIE DNA }

# COLOCAR ci, SOLVATACAO E SALVAR
saveamberparm CIED cie.top cie.crd
savepdb CIED cied.pdb
addions CIED Na+ 0
saveamberparm CIED cie_ci.top cie_ci.crd
solvatebox CIED TIP3PBOX 14
saveamberparm CIED cie_solv.top cie_solv.crd
savepdb CIED cie_solv.pdb

__________________

Foram realizadas três dinâmicas: model1, model4 e model5 (fazer). Os testes foram realizados no i7 e no cenapad

Modelo 4 (i7)

Em todos os testes o complexo se desligou do DNA. Motivos:

1) Erro na parametrização da carga do complexo(?)
2) DNA muito pequeno?
_____________________________
OBS:

Para rodar com o cuda, use primeiro o MPI na minimização e faça um pouco da dinâmica (md2.in) com MPI também. Assim quand tiver uma oa estrutura inicie o cálculo com gpu.


Observe que o complexo pode ficar emparelhado com o backbone com a região da piridina ou a região do oxindol. Observe que o =O o oxindol é bastante "pontiagudo" o que parece dificultar a interação com o DNA



Estudo 6DNA (12 bases) com a sequência AT ligado ao [Cu(isaepy(H2O)]+
Trajetória ~ ns
BOX: 24
(Cenapad - CUDA)
O sistema DNA apresentou bastante vibração. Ocorreu a separação da última sequência AT. A Timina se ligou om a adenina através do =O que não faz HB com o NH2 da adenina. Ocorreu também a interação entre bases de planos distintos.
O complexo foi "docado" com o cobre para o lado do solvente. O complexo se manteve estável a menos de algumas vibrações com a tentativa de trocar de sequência. Não houve a possibilidade de sair do DNA.


Estudo 3DNA (6 bases) com a sequência GC ligad ao [Cu(isaenim)]+ (forma enol).
Trajetória de ~56 ns
(Cenapad - CUDA)
BOX: 27
O composto inicialmente foi 'docado' com o cobre voltado para o DNA. O composto se desliga ficando preso pelo amino do imidazol. Em seguida o composto se vira e se liga com o lado hidrofóbico.

O DNA se mantém bem estruturado, sem grandes variações estruturais.


Estudo 6DNA (12 bases) com a sequência CG ligado ao [Cu(isaepy(H2O)]+
Trajetória ~ 54 ns
BOX: 24

(Cenapad: CUDA)

O complexo fica muito tempo preso ao backbone através da 5 posição de coordenação do cobre. Neste composto deixei uma água sempre fixa ao complexo. A última base também se desligou temporariamente.

Estudo 3DNA (6 bases) com a sequência CG ligado ao [Cu(isaepy(H2O)]+
Trajetória ~  ns
BOX: 24

(Cenapad: CUDA)

O composto fica bem mais instável, se desligando e religando ao DNA diversas vezes

O composto tenta se inserir entre as bases (intercalação) através da piridina. Também ocorre a interação frequente com o backbone através da 5 posição do cobre

A interação é mais frequente com o lado hidrofóbico, lembrando que neste modelo há uma água fixa ao composto.

[Cu(isaepy)]+ cotovelos

Comparação entre os resultados do docking para diferentes 'cotovelos'

Simulação feita com o Autodock vina e este atribui o seu campo eletrostático sem considerar as cargas fornecidas pelo usuário.

Observe que ainda não fiz tais testes com a dinâmica. E as que fiz foi com a forma 1

Resultados do docking da Forma 1:

Resultados Docking Forma 2:

Em resumo 
a Forma 1 pode intercalar pelo oxindol do lado do cobre ou pela piridina do lado invertido

a Forma 2 pode intercalar pela piridina pelo lado do cobre ou pelo oxindol do lado invertido.




Utilizei as Gasteiger partial atomic charges e como não há cargas para o cobre adicionei carga +1. Abaixo segue a atribuição das cargas para os sítios mais reativos. Entretanto como pode ser lido abaixo, os metais são tratados como doadores de ligação de hidrogênio e estas cargas são ignoradas no cálculo vina.

Cargas obtidas através do autodock
1)

ATOM      3 CU   LIG d   1       4.600   8.650   2.925  0.00  0.00     1.000 Cu
ATOM      1  N   LIG d   1       6.507   8.273   2.448  0.00  0.00    -0.266 N 
ATOM      2  N   LIG d   1       4.527   7.352   4.429  0.00  0.00    -0.284 N 
ATOM     19  N   LIG d   1       1.429   8.700   5.283  0.00  0.00    -0.220 NA
ATOM     20  O   LIG d   1       2.874   9.254   3.523  0.00  0.00    -0.168 OA



2)

ATOM      1 CU   LIG d   1       4.487   8.832   3.005  0.00  0.00     1.000 Cu
ATOM      2  N   LIG d   1       5.870   7.908   1.890  0.00  0.00    -0.266 N 
ATOM      4  N   LIG d   1       4.425   7.427   4.411  0.00  0.00    -0.284 N 
ATOM      3  N   LIG d   1       1.614   9.025   5.712  0.00  0.00    -0.220 NA
ATOM      5  O   LIG d   1       2.954   9.599   3.877  0.00  0.00    -0.168 OA

"AutoDock Vina ignores the user-supplied partial charges. It has its own way of dealing with the electrostatic interactions through the hydrophobic and the hydrogen bonding terms."

"Following Xscore, we formally treat metals as hydrogen bond donors"