From: Nassar, Omneya (omnassar_at_UTMB.EDU)
Date: Thu Sep 13 2018 - 21:54:08 CDT
Dear Ali,
Thank you for your reply and comments. I would like to first apologize for all the mistakes I made and for the poor formatting of this email. Here are my responses to your comments:
1. You’re absolutely right. By dipole I didn’t mean dipole moment and I meant two oppositely charged point charges (an electric dipole and not a permanent dipole).
Response: I’m sorry but I have never heard of an electric dipole and permanent dipole. There is only a dipole.
1. I agree with you that using tcl makes more sense but I was simply following NAMD tutorial on FEP calculations and they used pgn extension in their alchemical free energy calculation tutorial (Charging a spherical ion).
Response: You are correct. I pulled out my FEP tutorial and did see that they used “setup.pgn”. I apologize for the confusion.
1. So the setup.pgn (setup.tcl) script that they use in the tutorial for charging ONE atom is:
package require psfgen
topology ../common/top_all22_prot.inp
segment SOD {
residue 1 SOD
}
writepsf setup.psf
writepdb setup.pdb
Response: You are correct. I checked the tutorial and did see that they used this code. My apologies.
1. which generates setup.pdb:
REMARK original generated coordinate pdb file
ATOM 1 SOD SOD 1 q 0.000 0.000 0.000 -1.00 0.00 SOD NA
END
Response: This is almost correct. According to the fep tutorial, there is no “q”, no “NA” and -1 should be 1. The second line should look like this:
ATOM 1 SOD SOD 1 0.000 0.000 0.000 1.00 0.00 SOD
I apologize, but I previously confused the occupancy column with the beta column. The beta column comes after the occupancy column and before the segment column. Thus, having a -1 value in the occupancy column does not make sense because you must have the atom be present at the provided xyz coordinates a fraction number of times. For example, if you have “1” in the occupancy column, it means that the atom is located at that xyz coordinate for 100% of the time. If you have “0.5” in the occupancy column, it means that the atom is present at that specific xyz coordinate for 50% of the time. Thus, having “-1” in the occupancy column doesn’t make sense because how can an atom be present at a specific coordinate for -100% of the time? The psfgen program typically puts 1.00 in the occupancy column because that is the default coordinate assigned to the atom and only from structures that came from x-ray crystallography, NMR, Cryo-EM, etc have occupancy values that may deviate from 1.
I would like to clarify the composition of the pdb file. According to the NAMD tutorial, the pdb file columns “… in order from left to right are the record type, atom ID, atom name, residue name, [chain ID], residue ID, x, y, and z coordinates, occupancy, temperature factor (called beta), segment name, and line number." However, if the pdb file was “loaded into VMD and then written out as a new file, most of the extra information will be removed, the HETATM records will become ATOM records, and the previously empty chain ID field (between residue name and residue ID) will be set to X (unless present in the original file), and the line number will be omitted. ”
5. ATOM=ATOM ** atom number = 1 ** atom name = SOD ** residue name = SOD ** residue number = 1 ** chain ID = ??(NA?) ** xyz coordinates = 0.000 0.000 0.000 ** beta = -1 ** occupancy = 0.00 ** segment name = SOD ** NA?
Response: You are correct, except the occupancy = 1, beta = 0, segment name = SOD and periodic atom name = Na (not NA). I explain Na in response #10. I am not sure where the q and -1 came from.
6. and setup.psf file:
PSF
3 !NTITLE
REMARKS original generated structure x-plor psf file
REMARKS topology ../../common/top_all22_prot.inp
REMARKS segment SOD { first NONE; last NONE; auto angles dihedrals }
1 !NATOM
1 SOD 1 SOD SOD SOD 1.000000 22.9898 0
0 !NBOND: bonds
0 !NTHETA: angles
0 !NPHI: dihedrals
0 !NIMPHI: impropers
0 !NDON: donors
Response: I apologize, but I made a mistake when explaining the columns of the !NATOM section. According to the NAMD tutorial, the columns are "atom ID, segment name, residue ID, residue name, atom name, atom type, charge, mass and an unused 0." Thus, there is no section for the radius of the atom. This psf file looks correct.
8. as you can see, in the setup.pdb file the chain ID field is missing (or has been shifted to the end which doesn’t make sense!).
so the way I created the setup.pdb and setup.psf files for my simulation was to modify setup.pgn to have two sodium atoms instead of one:
package require psfgen
topology ../common/top_all22_prot.inp
segment SOD1 {
residue 1 SOD
}
segment SOD2 {
residue 1 SOD
}
writepsf setup.psf
writepdb setup.pdb
Response: I see your logic in setup.pgn and it makes sense. I apologize for causing any confusion. The way I usually build pdb files with multiple segments is by having a pdb file for each segment (e.g. segment SOD1 {pdb file1.pdb} // segment SOD2 {pdb file2.pdb}). I found that this method makes it easier to debug any errors that may arise (even if the psfgen plugin runs smoothly) because I can trace back every character of the new pdb file back to the old pdb file.
9. by doing so and using the command source setup.pgn in VMD TkConsole, setup2.pdb and setup2.psf files are generated as follow:
Response: It looks like you made a typo when writing out the file names. Instead of "setup2.pdb" and "setup2.psf", it should be "setup.pdb" and "setup.psf", respectively. Other than that, using the command source setup.pgn is absolutely correct.
10. setup.pdb:
ATOM 1 SOD SOD 1 0.000 0.000 0.000 -1.00 0.00 SOD1 NA
ATOM 2 SOD SOD 1 0.000 0.000 0.000 -1.00 0.00 SOD2 NA
Response: Once again, you cannot have -1.00 in the occupancy column. VMD puts a "1.00" for the occupancy, not -1.00. I have had instances when the periodic atom name (e.g. C for carbon, N for nitrogen, O for oxygen, H for hydrogen, etc.) is located in the same location as the NA in my previous pdb files, so that may be where the NA is coming from. However, it should be "Na" instead of "NA", right? From my experience, VMD confused one of the phosphate atoms as "Pb" (lead) instead of "P" (phosphate), and I found this mistake in my pdb file on the same column you have "NA" in. However, note that they use the standard nomenclature and capitalization for the periodic atoms (the first letter is capitalized and the second letter is lowercase).
11. setup.psf
PSF
4 !NTITLE
REMARKS original generated structure x-plor psf file
REMARKS topology ../common/top_all22_prot.inp
REMARKS segment SOD1 { first NONE; last NONE; auto angles dihedrals }
REMARKS segment SOD2 { first NONE; last NONE; auto angles dihedrals }
2 !NATOM
1 SOD1 1 SOD SOD SOD 1.000000 22.9898 0
2 SOD2 1 SOD SOD SOD 1.000000 22.9898 0
0 !NBOND: bonds
0 !NTHETA: angles
0 !NPHI: dihedrals
0 !NIMPHI: impropers
Then I modified the setup.pdb file and change the x coordinate of the two atoms to (-r/2) and (r/2) where are is the desired distance between the two atoms (apologies for using d instead of r). Also modifying the setup.psf file by changing the charge of the sodium atom in segment S1 from 1.00000 to -1.00000.
Response: Everything you said here is correct except for the typo "S1" instead of "SOD1".
12. package require solvate
solvate setup.psf setup.pdb -minmax {{-20-(r/2) -20-(r/2) -20-(r/2)} {20+(r/2) 20+(r/2) 20+(r/2)}} -o solvated
Response: You have to specify the value of r in the tcl script. Otherwise, the code will read it as an undefined variable and not know what to do with it. Other than that, everything else here looks correct. Optional: if you want to have a well solvated box that prevents periodic images of the sodium ions from interacting with the sodium ions, a 40 x 40 x 40 angstrom^3 solvated box will be sufficient enough to do so and you will not have to include r/2. However, having so many waters in the box may be excessive and more computationally expensive than necessary.
13. which generates solvated.pdb (for r=3A):
Response: This looks correct
14. and lastly changing the beta column for the two SOD atoms from 0.00 to 1.00 and finally saving as solvated.fep.
Note that for some reasons the residue names for the two SOD atoms have changed to S (instead of SOD) which apparently doesn’t cause a problem either!
Response: S and W look like they are chain IDs, not residue names. The residue names are still there (SOD and TIP3).
15. "The fields in the atom section are atom ID, segment name, residue ID, residue name, atom name, atom type, charge, mass, and an unused 0.”
Response: This is correct and I mentioned this above. I apologize for not reading this earlier, but I read and responded to each of your comments by each section. Thank you for correcting me and I apologize for saying something wrong.
16. I checked the solvate.pdb file in NAMD FEP tutorial and there’s no spacing between TIP3 and W (however I agree that it should be) and that doesn’t cause a problem in their simulation and in mine as well (my simulations don’t crash).
Response: You are correct and I am sorry for the confusion.
17. You’re absolutely right. All beta values of Water molecules are automatically set to zero in the fee file.
Response: You mean the fep file, not fee file?
18. Yes! I meant the effects of the interaction of the sodium ion with it’s periodic image. I understand the Hummer et al. correction for the situation that we have one ion in the simulation box. But in my case that the entire system presumably remains neutral the whole time, I was wondering what would be the correction?
All my simulations ran with no problem but I need to know what’s the right correction term term to add to the raw energies obtained from FEP.
Response: You should not have the sodium ions interacting with their periodic image; that's why many people solvate simulation boxes with 10 angstroms or more padding of water. You are trying to calculate the electrostatic free energy of charging a sodium ion to +1 and charging a sodium ion to -1. You are not calculating the interaction of a sodium ion with it's periodic image; that would be like calculating the interaction of two ions with the same charge. Also, as I mentioned previously, you do not need to include a correction to the free energy because you system remains neutral throughout the entire simulation. Thus, the particle mesh ewald method will be able to handle the long range electrostatics properly. Finally, total free energy of your simulation will be zero because the electrostatic free energy of one sodium ion will be equal and opposite in sign to the electrostatic free energy of the other sodium ion. These values will cancel and you will get zero.
I hope my responses were helpful and I apologize for any mistakes or redundancies I made in any of my responses.
Thank you,
Omneya Nassar
PhD Candidate in Pharmacology and Toxicology
University of Texas Medical Branch
301 University Blvd
Galveston, TX 77555
E-mail: omnassar_at_utmb.edu
Phone: 409.772.0731
________________________________
From: Ali Mehdizadeh Rahimi [mehdizadehrahimi.a_at_husky.neu.edu]
Sent: Thursday, September 13, 2018 10:58 AM
To: Nassar, Omneya
Subject: Re: namd-l: Ewald size-dependence correction in FEP simulation of an electric dipole
WARNING: This email originated from outside of UTMB's email system. Do not click links or open attachments unless you recognize the sender and know the content is safe.
Dear Omneya,
Thank you so much for taking the time and posting a detailed answer to my question. I would like to further discuss some issues with you and resolve some confusions that I have:
1. A dipole moment can only happen between two atoms covalently bonded and the difference between the partial charges of each atom is 0.5 or greater. Because sodium ions form ionic bonds instead of covalent bonds, their interactions are purely electrostatic and van der Waals interactions (nonbonded)--_000_EE608134F0CB1346B16AC09716744C16948DACF6GRMBX1utmbedu_--
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2019 - 23:20:12 CST