Research:amber prog tips

VALENCE

  opbend is included in angle energy

 in amoeba_interface.f:

  call AM_BONDS_eval(crd,frc,ene(1),vir_tensor)
  call AM_UREYB_eval(crd,frc,ene(2),vir_tensor)
  call AM_REG_ANGLES_eval(crd,frc,ene(3),vir_tensor)
  call AM_TRIG_ANGLES_eval(crd,frc,ene(4),vir_tensor)
  call AM_OPBEND_ANGLES_eval(crd,frc,ene(5),vir_tensor)
  call AM_TORSIONS_eval(crd,frc,ene(6),vir_tensor)
  call AM_PITORSIONS_eval(crd,frc,ene(7),vir_tensor)
  call AM_STRETCH_BEND_eval(crd,frc,ene(8),vir_tensor)
  call AM_TOR_TOR_eval(crd,frc,ene(9),vir_tensor)
  call AM_STRETCH_TORSIONS_eval(crd,frc,ene(10),vir_tensor)

AND

ebond = ene(1)!bonds energy
eangle = ene(2)+ene(3)+ene(4)+ene(5)+ene(8)!angle energy
etors = ene(6)+ene(7)+ene(9)+ene(10)!torsion energy

VDW short range ADJUST & softcore

(a test dire: /home/pren/AMOEBA/solvation/cresol/vdw0.4)

For vdw, the AM_VDW_DIRECT_ene_frc_i subroutine, where the softcore mod is included, compute interaction for non-adjusted pair.

So it is important to understand how the adjust list is created.

In *.prmtop, the AMOEBA_ADJUST_LIST gives all the pairs in the adjust list in (3xNUM) arrays:

 1 2 5 means atom 1, 2 pair and adjust type 5.

The waight for type 5 is defined in AMOEBA_ADJUST_VDW_WEIGHTS_LIST of prmtop. 5 is 1-2, 6 is 1-3, 7 is 1-4...

I have looked at the amoeba_parm.f amd amoeba_vdw.f carefully, it turns out the assumption that the soft-core mod only applies to intermolecular pairs is not correct. The AM_VDW_DIRECT_ene_frc_i routine go through any atom pairs that are not in the adjust list – this where the softcore modification. There is another routin in amoeba_vdw.f AM_VDW_ADJUST_ene_frc, this where we assumed that intramolecular interaction goes so that no soft-core needed here.

In fact the amoeba_parm.f shows that the adjust list (or exclusion list) is built from both polarization group and 1-2 (bonded), 1-3(angle)…1-5 interactions. So any atom belongs to these categories, will be in the adjust (exclusion) list. For benzamidines, cresol, isoprop, the whole molecule is a polarization group so any vdw pair is in this list. So the current softcore code is correct.

But for NMA (3 polarization groups), or any other large molecule that has more than one polarization group, the soft-core will affect intramolecular pairs as well. For NMA, the pairs are between H aotms in the two CH3 groups, which belong to different pol groups and are 1-6 pairs. If you run analyze on a single NMA molecule, you can calculate these contributions which are 0.06 kcal/mol total. So for NMA the code is wrong but didn’t have much effect. Oscar, can you confirm this contribution is indeed small? In /home/pren/AMOEBA/solvation/cresol/vdw0.4 I set up an example for a single NMA, and print out the i j pairs from amoeba_vdw.f. You can clearly see the atom pairs (the top half is printed from the soft-core loop; the atom id matches nma.xyz). we can certainly grow this intra-mol contribution back in gas-phase if needed.

local frame

In /home/pren/DARDEN_amber/SRC/build_amoeba/amoeba_parm.f90 the tinker multipoles of z-then-x and bisector are recorded as frame defining vectors in subroutine "AM_PARM_get_multipole"

2013 APril

4/13/2013 in tinker empole.f 2009: if (i14(j,ii) .eq. ip11(k,ii))

    &            pscale(i14(j,ii)) = 0.5d0 * pscale(i14(j,ii)) 

after 2010-sep:

if (i14(j,ii) .eq. ip11(k,ii))
    &            pscale(i14(j,ii)) = p4scale * p41scale

With p41scale default to 1.0

In build_amoeba/amoeba_parm.f:

! note polar wt of 1/2 for pair in polar group and 1-4

 polar_wt(7) = 0.5d0; polar_wt(8) = 1.d0; polar_wt(9) = 1.d0

---

4/2013 pmemd-10: softcore already changed to the soft_atm.txt

2013 The ligand_idx1 is obsolte. By default pmemd look for a file soft_atm.txt Note soft_lambda replaces soft_lamda http://water.bme.utexas.edu/wiki/index.php/Software:amber#Instruction

Distance restraint MROPT no longer working?

gvim -d  ~/pmemd-oscar/src/nmr_calls.f90 ~/DARDEN_amber/pmemd/src/

Bnd constraint no longer work in the nmr control of sander or pmemd.

sander reprocess trajectory

see amoeba-peptide/fep-test
imin=5 will call trajene and minimize/runmin.f (maxcycle=1) every frame
Changed dynlib.f 10f8.3 to 8f10.4, same with trajene, to 8f10.4
./sander -i pepc.min -c pepc.inpcrd -x mdcrd -p pepc.prmtop -o pepc.minout -r pepc.restrt
-x mdcrd specify trajctory generated previously, change the pepc.prmtop rerun sander to evaluate energy at perturbed parameters

---

ntwx=1 in your input Try to set ntwx=1 in your input file. The default value is 0 and then the trajectory file is not opened.

you didn't specify the trajectory using the -x flag.
-c and -ref should still be input coordinates, not
trajectory. -ref is only needed if you are using positional
or RMSD restraints.
look at tests/trajene for an example.

>I am trying to run sunder with imin=5 parameter for trajectory analysis. But >sander cannot process my trajectory file. Sander reads rst or inpcrd files but >gives out this message when I try to run with a mdcrd file: > >fmt: read unexpected character >apparent state: internal I/O >last format: (i6,e15.7) >lately reading sequential formatted internal IO >./min5_EN.q: line 7: 15353 Aborted $AMBERHOME/exe/sander -O -i >min5_EN.in -o min5_EN.out -p 3_20.prmtop >-c 3_20.mdcrd -r new.rst -ref 3_20.mdcrd > >I realized that there is a number at the beginning of my inpcrd files (something >like 39150) and that does not exist in a mdcrd file. I suspect that may be a >reason. > >Is there an intermediate step before analyzing a trajectory file with imin=5? > >also my sander input file is: >minimization_step1 > &cntrl > imin = 5, > maxcyc = 1, > ntb = 1, > ntr = 0, > cut = 10 > / >END >

---

Convert mdcrd

convert mdcrd (can be viewed by vmd directly) or incrd to pdb

/opt/amber7/exe/ptraj pepc.prmtop pepc.scrpt

pepc.scrpt:
       trajin mdcrd 
       trajout tra.pdb pdb 

for octahedron:
trajin mdcrd 100 150 10
trajout oct.pdb pdb
image origin familiar

convert to pdb for amber simulation

There seems to be problem with your ca/ef hand run, the atom ordering in xyz (and .dyn) and pdb are not the same. D id you try diff to check?

For example look at atom 208 in both files.

This is due to xyzpdb (siwth seq present) changed the order within each residue. We probably have to do some manual

editing.

Here is what that will work - in the future we should follow: 1) starting with pdb. Pdbxyz create xyz with parameters, and stripped ligand/water. A seq file is created at the sa me time 2) add ligand soak water, run minimization and MD with TINKER. Dyn file is created 3) do xyzpdb on the box.xyz (cp from last dynamics frame) with seq. The atom order will be changed from xyz file w ithin each residue, but everything else is fine 4) repeat but without seq file. The order will be fine but many others are wrong. 5) Take the protein atom name column from 4) and replace the column in 3) 6) notice by doing 5) the atom name order is changed to that of xyz but the coordinates are wrong. This is fine sin ce the real coordinates will be read from dyn file

Try this with your efhand see if this make sense.