Initial pdb file processing: histidines.inp

Before starting with the minimization and solvation, one ought to decide on the protonation states of the titratable groups. Often this means analyzing what is around each of these groups after carefully reading the pdb remarks to see if there are any hints about whether they are coordinated. If they are coordinated, then it is important to check the distances in the coordination complex to decide which groups should or should not be protonated.
Run charmm making sure to identify the pdbfile (1fsc.pdb):
% charmm pdbfile=1fsc < histidines.inp
Two methods for looking aqt what the protonation state of the histidines are presented here and both require the user to analyze the results and to make the final decision. This example has 2 histidines (resid 6 and 29).
COOR MINDist (input script) will give the closest atom to each ring nitrogen in each histidine (actual output file) important output as follows. Residue 6 has an electron donor near ND1 indicating it should be protonated, and NE2 has an electron acceptor near it indicating it should not be protonated. We should decide to make residue 6 a HSD. We will make residue 29 a HSE since there is nothing near ND1 but an acceptor near NE1.
            MINIMUM DISTANCE FOR ALL SELECTED ATOMS
   83 1FSC HSD  6    ND1 - 106 1FSC THR  7    O         2.8497
  
            MINIMUM DISTANCE FOR ALL SELECTED ATOMS
   88 1FSC HSD  6    NE2 - 583 1FSC ARG  37   HH12      1.9869
  
            MINIMUM DISTANCE FOR ALL SELECTED ATOMS
  446 1FSC HSD  29   ND1 - 420 1FSC ARG  28   HB1       3.6827
  
            MINIMUM DISTANCE FOR ALL SELECTED ATOMS
  451 1FSC HSD  29   NE2 - 690 1FSC ASN  47   OD1       2.6855
  
Histidine HSD resid 6
Histidine HSE resid 29
COOR DIST will give all atoms within a cutoff (here we use 3 A) and the followong command (input script ) prints a list of atoms near the histidine ring nitrogens (NE2 and ND1) (actual output file)
 
 CHARMM>    coor dist cut 3.0 sele resn HS* .and. (type NE2 .or. type ND1 ) end -
 CHARMM>                           sele .not. resn HS* end
 SELRPN>      4 atoms have been selected out of    906
 SELRPN>    872 atoms have been selected out of    906

            DISTANCES FOR SELECTED ATOMS
   83 1FSC HSD  6    ND1 - 106 1FSC THR  7    O         2.8497
   88 1FSC HSD  6    NE2 - 574 1FSC ARG  37   HG2       2.8153
   88 1FSC HSD  6    NE2 - 576 1FSC ARG  37   HD1       2.8169
   88 1FSC HSD  6    NE2 - 581 1FSC ARG  37   NH1       2.8944
   88 1FSC HSD  6    NE2 - 583 1FSC ARG  37   HH12      1.9869
  451 1FSC HSD  29   NE2 - 685 1FSC ASN  47   HA        2.9270
  451 1FSC HSD  29   NE2 - 690 1FSC ASN  47   OD1       2.6855
 TOTAL EXCLUSION COUNT =         0
 TOTAL 1-4 EXCLUSIONS  =         0
 TOTAL NON-EXCLUSIONS  =         7
  
 CHARMM>     
Since the decision has been made to change one HSD to a HSE, the appropriate patch needs to be added to the solvate script before the psf is finalized and written. The patch is:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Patch for changing HSD to HSE for this particular case
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
patch hs2  @pdbfile 29
rename resn HSE sele resi 29 end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      Done Patch
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!