Corso di Laurea Magistrale in Scienze Chimiche, Curriculum "Struttura, dinamica e reattività chimica"

Struttura e Dinamica Molecolare di Sistemi Biologici - 2010/2011

Laboratorio 3 - (6 dicembre 2010)

Simulazioni con solvente.

Preparazione input

  1. Volume cella e densità solvente

    1. lato cella: è opportuno che sia almeno il doppio della lunghezza della proteina
    2. densità dell'acqua:

      \begin{eqnarray*} \rho/amu\cdot\AA^{-3} & = & \rho/g\cdot cm^{-3}\cdot\frac{g}{a... ...)^{3}\\ & = & \rho/g\cdot cm^{-3}\cdot0.6023\\ & \simeq & 0.6\end{eqnarray*}

      Per $m^{3}$ molecole di acqua in un box di lato $L$ Angstrom:

      \begin{displaymath} \rho/amu\cdot\AA^{-3}=0.6=\frac{18m^{3}}{L^{3}}\end{displaymath}

      da cui

      \begin{displaymath} L=3.107m\end{displaymath}

      Ad esempio

      \begin{eqnarray*} m & = & 10\\ L & = & 31.07\end{eqnarray*}

    3. Quando si inserisce il soluto, questo si sovrappone al solvente; si cancellano molecole di solvente entro un certo raggio dagli atomi del soluto, raggio definito dal parametro INSERT del blocco &SOLVENT (due molecole si considerano sovrapposte se la somma delle loro distanze è minore del loro raggio di Lennard-Jones moltiplicato per INSERT). Notare che in questo modo la densità (e pressione) sperimentale non è più rispettata: la soluzione ideale sarà fare simulazione a pressione costante, almeno all'inizio, e lasciare che il volume si aggiusti (cfr. avanti).
  2. creare sulla propria directory di lavoro il file di input, usando come modello ad es:

    ~signorini/biomol/orac/chigno-slv/slv-nvt.in

  3. copiare nel proprio ambiente di lavoro il file di coordinate della molecola di $H_{2}O$:

    ~signorini/biomol/orac/pdb/water.pdb

  4. Notare le differenze del file di input rispetto al file utilizzato per la simulazione nel vuoto (magari usando il comando ediff di Emacs ):

    1. aggiungere solvente
    &SETUP
    CRYSTAL 31.07 31.07 31.07
    &END

    &SOLUTE
    COORDINATES chigno-01.PDB
    &END

    &SOLVENT
    CELL SC

    INSERT 0.8

    COORDINATES ../pdb/water.pdb

    GENERATE RANDOMIZE 10 10 10
    &END

    &PARAMETERS
    [...]

    JOIN SOLVENT
    spce
    END
    &END
    1. aggiornare i time-step
    &INTEGRATOR
    TIMESTEP 10.0

    MTS_RESPA
    step intra 2

    step intra 2

    step nonbond 2 4.7

    step nonbond 3 7.5 reciprocal

    step nonbond 1 9.7

    test_times OPEN energie
    END
    &END
    1. Ewald
    &POTENTIAL
    EWALD PME 0.43 30 30 30 4

    [...]
    &END
    1. salvare un restart
    &INOUT
    RESTART
    write 500.0 OPEN chigno.rst
    END

    ASCII 200.0 OPEN chigno.pdb

    [...]
    &END

Verifica conservazione E

È bene verificare la scelta dei timestep e dei cutoff associati. La quantità che si deve conservare è ETOT nel file energie. Notare che se il solvente è stato appena creato l'eccesso di energia è grande e può darsi che il termostato non ce la faccia a assorbirlo: bisogna prima fare un'equilibratura di qualche picosecondo.

Equilibrazione

  1. Anche in simulazioni a $T=cost$ conviene fare una prima parte di equilibrazione, per evitare che il termostato si scaldi troppo. Eventualmente si riparte dal restart.

    Figura: confronto di energia cinetica, totale e potenziale con e senza ciclo di equilibratura

    \includegraphics{nvt-rej+norej}

  2. notare quanto è più lenta la simulazione aggiungendo il solvente (circa 10 volte)
  3. monitorare le energie, in particolare:

    1. Energia totale EREAL [scende , poi cost]
    2. EPTOT
    3. EKIN [oscilla intorno a 300K]
    4. ESLV [scende]
    5. ESLV-SLT [scende?]
    6. ESLT [potrebbe salire]
    Tutte queste grandezze sono listate nel file energie.

    Si controllano usando gnuplot

  4. Visualizzare il sistema in VMD

    1. vedere se la RMSD rispetto alla conformazione sperimentale è minore che nel vuoto

usare pressione costante

È opportuno fare una simulazione $N,p,T$ per avere la densità corretta. Si parte da un alto valore del raggio di sovrapposizione e si permette alla cella di comprimersi sotto l'effetto della pressione esterna ($0.1MPa=1atm$)

&SOLVENT
CELL SC

INSERT 1.4

COORDINATES ../pdb/water.pdb

GENERATE RANDOMIZE 10 10 10
&END

&SIMULATION
...

THERMOS
...
END

ISOSTRESS PRESS-EXT 0.1 BARO-MASS 10.
&END

Figura: confronto di energia cinetica, totale e potenziale in NVT e NpT

\includegraphics{npt+nvt}

In questa simulazione monitorare le energie delle variabili estese:

Simulazioni di solvente puro

\usebox{\@tempboxa}