CHARMM Development Project Forums User Discussion & Questions Energy terms, Constraints, Restraints, and Solvation Can Charmm compute pi stacking interaction?

For a system of protein containing phenylalanine residue, I used following command to calculate pi-pi interaction
SKIP ALL EXCL vdw Elec
interaction select resname phe .and. type C* end sele resname phe .and. type C* end
ECHO "PHEC PHEC" @X ?VDW ?ELEC

I expected to see at least one of the two terms to be negative. But both VDW and ELEC are positive.

You select only the carbon atoms, so the electrostatic interaction will come out as repulsive since all the ring-carbons have the same charge. Try with a selection that includes all the ring atoms, and no backbone atoms. Something like: sele resname phe .and. (chem ca .or. chem hp) end. You may also want to look at specific phe-phe pairs. If you want to play with it, you can easily create a system with just two phenylalanine rings and move them apart while evaluating the energies every 0.1Ã… or so; perhaps also with the rings in different orientations

I selected all the carbons and hydrogens of ring atoms from PHE residue as suggested. Still, I got positive VDW and ELEC term. Since charges of CA and HP are opposite, should I simply compute the interaction between CA and HP atoms instead ?(ex: interaction select resname phe .and. chem CA end sele resname phe .and. chem HP end).

Afterwards, I did get negative VDW and ELEC by computing interaction between CA and HP. And I also computed interaction between CA and CA and between HP and HP. But the sum of positive interaction is larger than the negative CA -HP. This result still does not make sense to me since in reality the Pi stacking interaction takes account of every atom of PHE ring.

From VMD, the PHE rings of crystal should be close enough.

OK. The relevant aspect of the energies computed by CHARMM is that they should always be used in a relative sense, so you should always compute the difference wrt some reference state. This means that even though the energy values you compute are all positive, there may still be an attraction between the two rings, as long as there is a minimum in the energy surface. If I remember correctly the minimum conformation for two benzenes using the CHARMM energy function is where the edge of one interacts with the flat surface of the other (more or less).

There are also interactions within each ring, which are probably positive, but not important for the interaction between the rings. To avoid this complication you can use the interaction_energy command (energy.info).