added LJ computations

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.1-DEV"
 [deps]
 CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 Chemfiles = "46823bd8-5fb3-5f92-9aa0-96921f3dd015"
+FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MolecularMinimumDistances = "46954fcc-1dd7-4352-b01d-a39240e91dc8"
 PDBTools = "e29189f1-7114-4dbd-93d0-c5673a921a58"

README.md

@@ -25,35 +25,37 @@ using SolventShellInteractions
 dir="./test/files"
 
 # read pdb file
-pdb = readPDB("$dir/simulacao_EMIMDCA.pdb")
+pdb = readPDB("$dir/system.pdb")
 
 # solute atoms
 solute = select(pdb, "protein")
 solvent = select(pdb, "resname EMI")
 
 # trajectory file (Gromacs xtc only)
-trajectory = "$dir/simulacao_EMIMDCA_curta.xtc"
+trajectory = "$dir/simulation_short.xtc"
 
-# topology files
-top_files = [ "$dir/topol.top", "$dir/tip3p.itp" ]
+# topology file
+topology_file = "$dir/processed.top"
 
 # distance of the first dip in the distribution
 cutoff = 10.
 
 # compute electrostatic potential
-u = electrostatic_potential(
+q, lj = nonbonded(
     solute,
     solvent,
     cutoff,
     trajectory, 
-    top_files,
+    topology_file,
 )
 
 plot(
-    u,
+    [ q lj ],
+    labels=[ "electrostatic" "Lennard-Jones" ],
     xlabel="step",
-    ylabel="electrostaic potential / kJ / mol",
-    linewidth=2, framestyle=:box, label=nothing
+    ylabel="potential energy / kJ / mol",
+    linewidth=2, 
+    framestyle=:box,
 )
 ```
 
@@ -69,6 +71,7 @@ Three keyword arguments can be passed to the `electrostatic_potential` function:
 | keyword |  values | default |  meaning  | 
 |:-------------:|:---------------:|:--:|:-------------|
 | `show_progress` | `true/false` | `true` | Display or not the progress bar.  | 
+| `combination_rule` | `:geometric/:arithmetic` | `:geometric` | Type of combination rule for the sigma parameters of the LJ potential.  | 
 | `standard_cutoff` | `true/false` | `false` |Use a standard cutoff, and thus no interaction above the cutoff distance will be considered. | 
 | `shift` | `true/false` | `false` | Use a shifting function, as described [here](https://www.ks.uiuc.edu/Research/namd/2.10/ug/node23.html).  | 
 

src/SolventShellInteractions.jl

@@ -4,32 +4,36 @@ using LinearAlgebra: norm
 using MolecularMinimumDistances
 using StaticArrays
 using ProgressMeter
+using FastPow
 
 import CellListMap: wrap_relative_to
 import Chemfiles
 import PDBTools
 
 export FFType
-export read_ffcharges, read_ffcharges!
-export assign_charges
-export electrostatic_potential
+export read_forcefield, read_forcefield!
+export assign_forcefield
+export nonbonded
 
 Base.@kwdef struct FFType
     type::String
     name::String
     resname::String
     charge::Float64
+    eps::Float64
+    sig::Float64
 end
 
 # Read charges from topology file, for the first time
-function read_ffcharges(topology::String)
-    ffcharges = FFType[]
-    read_ffcharges!(ffcharges, topology)
-    return ffcharges
+function read_forcefield(topology::String)
+    ff = FFType[]
+    read_forcefield!(ff, topology)
+    return ff
 end
 
 # Append new data to the ffcharges array, given a new topology file
-function read_ffcharges!(ffcharges::Vector{FFType}, topology::String)
+function read_forcefield!(ff::Vector{FFType}, topology::String)
+    # Read charges and atom types
     reading = false
     open(topology) do file
         for line in eachline(file)
@@ -56,51 +60,100 @@ function read_ffcharges!(ffcharges::Vector{FFType}, topology::String)
                 type = type_data[2], 
                 name = type_data[5], 
                 resname = type_data[4], 
-                charge = parse(Float64, type_data[7])
+                charge = parse(Float64, type_data[7]),
+                eps = -1.,
+                sig = -1.,
             )
             # if the atom type is not repeated, add to the list
-            if isnothing(findfirst(isequal(fftype), ffcharges))
-                push!(ffcharges, fftype)
+            if isnothing(findfirst(isequal(fftype), ff))
+                push!(ff, fftype)
             end
         end
     end
-    return ffcharges
+
+    # Read LJ parameters, given classes
+    reading = false
+    open(topology) do file
+        for line in eachline(file)
+            # start reading atoms section
+            if occursin("[ atomtypes ]", line)
+                reading = true
+                continue
+            end
+            if !reading 
+                continue
+            end
+            # end of atoms section 
+            type_data = split(strip(line))
+            if length(type_data) == 0
+                reading = false
+                continue
+            end
+            # comment line on atom section
+            if type_data[1] == ";"
+                continue
+            end
+            ilast = findfirst(isequal(";"), type_data)
+            if isnothing(ilast)
+                ilast = length(type_data) + 1
+            end
+            for (itype, atom_type) in pairs(ff)
+                if atom_type.sig == -1
+                    if atom_type.type == type_data[1]
+                        ff[itype] = FFType(
+                            type = ff[itype].type, 
+                            name = ff[itype].name, 
+                            resname = ff[itype].resname, 
+                            charge = ff[itype].charge,
+                            eps = parse(Float64,type_data[ilast-1]),
+                            sig = 10*parse(Float64,type_data[ilast-2]),
+                        )
+                    end
+                end
+            end
+        end
+    end
+    return ff
 end
 
-function assign_charges(atoms, ffcharges; warn_aliasing = true)
-    charges = zeros(length(atoms))
+function assign_forcefield(atoms, ff; warn_aliasing = true)
+    params = Vector{typeof((q=0.,eps=0.,sig=0.))}(undef, length(atoms))
     for (iat, at) in pairs(atoms)
-        itype = findfirst( 
-            type -> ( type.resname == at.resname && type.name == at.name ),
-            ffcharges
-        )
+        itype = findfirst( type -> ( type.resname == at.resname && type.name == at.name ), ff )
         if isnothing(itype)
             # checking for name aliasing for H atoms
             if length(at.name) == 4
                 aliased_name = at.name[2:4]*at.name[1]
                 itype = findfirst( 
                     type -> ( type.resname == at.resname && type.name == aliased_name ),
-                    ffcharges
+                    ff
                 )
                 if !isnothing(itype)
                     warn_aliasing && println(" Warning: aliasing $(at.resname): $(at.name) to $(aliased_name).")
                 end
             end
             if isnothing(itype)
-                error("Could not find charges for atom: $(at.name) $(at.resname)")
+                error("Could not find parameters for atom: $(at.name) $(at.resname)")
+            else
+                if ff[itype].charge == -1
+                    error("Could not find charges for atom: $(at.name) $(at.resname)")
+                end
+                if ff[itype].eps == -1 || ff[itype].sig == -1
+                    error("Could not find Lennard-Jones parameters for atom: $(at.name) $(at.resname)")
+                end
             end
         end
-        charges[iat] = ffcharges[itype].charge
+        params[iat] = (q=ff[itype].charge, eps=ff[itype].eps, sig=ff[itype].sig)
     end
-    return charges
+    return params
 end
 
-function electrostatic_potential(
+function nonbonded(
     solute::AbstractVector{PDBTools.Atom}, 
     solvent::AbstractVector{PDBTools.Atom}, 
     cutoff::Real,
     trajectory::String,
-    topology_files::Vector{String};
+    topology_file::String;
     standard_cutoff::Bool = false,
     shift::Bool = false,
     show_progress::Bool = true,
@@ -114,32 +167,30 @@ function electrostatic_potential(
     natoms_per_molecule = length(PDBTools.select(solvent, "resnum $(solvent[1].resnum)"))
 
     # Read charges from topology files
-    ffcharges = FFType[]
-    for top_file in topology_files
-        read_ffcharges!(ffcharges, top_file)
-    end
+    ff = read_forcefield(topology_file)
 
     # Assign charges
-    solute_charges = assign_charges(solute, ffcharges; warn_aliasing = false)
-    solvent_charges = assign_charges(solvent, ffcharges; warn_aliasing = false)
+    solute_params = assign_forcefield(solute, ff, warn_aliasing = false)
+    solvent_params = assign_forcefield(solvent, ff, warn_aliasing = false)
 
     # Open trajectory with Chemfiles
     traj = Chemfiles.Trajectory(trajectory)
 
-    u = zeros(length(traj))
+    electrostatic_potential = zeros(length(traj))
+    lennard_jones = zeros(length(traj))
     show_progress && (p = Progress(length(traj)))
     for iframe in 1:length(traj)
         frame = read(traj)
         coor = reinterpret(reshape, SVector{3,Float64}, Chemfiles.positions(frame))
         unit_cell = Chemfiles.lengths(Chemfiles.UnitCell(frame))
 
         # Compute electrostatic potential
-        u[iframe] = electrostatic_potential(
+        electrostatic_potential[iframe], lennard_jones[iframe] = nonbonded(
             solute_indexes,
             solvent_indexes,
             natoms_per_molecule,
-            solute_charges,
-            solvent_charges,
+            solute_params,
+            solvent_params,
             coor,
             unit_cell,
             cutoff;
@@ -150,20 +201,31 @@ function electrostatic_potential(
         show_progress && next!(p)
     end
 
-    return u
+    return electrostatic_potential, lennard_jones
 end
 
-function electrostatic_potential(
+function lj(eps_x, eps_y, sig_x, sig_y, d; combination_rule = :geometric) 
+    if combination_rule == :geometric
+        sig = sqrt(sig_x * sig_y)
+    else
+        sig = 0.5 * (sig_x + sig_y)
+    end
+    lj = @fastpow 4 * sqrt(eps_x * eps_y) * ( (sig / d)^12 - (sig / d)^6 )
+    return lj
+end
+
+function nonbonded(
     solute::AbstractVector{Int}, 
     solvent::AbstractVector{Int}, 
     natoms_per_molecule::Int, 
-    solute_charges::AbstractVector{<:Real}, 
-    solvent_charges::AbstractVector{<:Real}, 
+    solute_params::AbstractVector{<:NamedTuple}, 
+    solvent_params::AbstractVector{<:NamedTuple}, 
     coor::AbstractVector{<:SVector}, 
     unit_cell, 
     cutoff::Real;
     standard_cutoff::Bool = false,
-    shift::Bool = false
+    shift::Bool = false,
+    combination_rule = :geometric
 )
 
     # Define simulation box in this frame
@@ -182,6 +244,7 @@ function electrostatic_potential(
     # molecules that have some atom within the cutoff 
     nbatches = Threads.nthreads()
     electrostatic_potential = zeros(nbatches)
+    lennard_jones = zeros(nbatches)
     Threads.@threads for ibatch = 1:nbatches
         for isolvent = ibatch:nbatches:length(list)
             # skipt if the solvent molecule is not within the cuotff
@@ -192,23 +255,26 @@ function electrostatic_potential(
             ilast = ifirst + natoms_per_molecule - 1
             for iatom_solvent = ifirst:ilast
                 x = xsolvent[iatom_solvent]
-                qsolvent = solvent_charges[iatom_solvent]
+                qx, eps_x, sig_x = solvent_params[iatom_solvent]
                 for (iatom_solute, y) in pairs(xsolute)
-                    qsolute = solute_charges[iatom_solute]
+                    qy, eps_y, sig_y = solute_params[iatom_solute]
                     y = wrap_relative_to(y,x,box) 
                     d = norm(y-x)
                     if d < box.cutoff || !standard_cutoff
-                        qpair = qsolvent * qsolute / d
+                        qpair = qx * qy / d
+                        ljpair = lj(eps_x, eps_y, sig_x, sig_y, d; combination_rule = combination_rule)
                         if shift
-                            qpair -= qsolvent * qsolute / box.cutoff
+                            qpair -= qx * qy / box.cutoff
+                            ljpair -= lj(eps_x, eps_y, sig_x, sig_y, box.cutoff; combination_rule = combination_rule)
                         end
                         electrostatic_potential[ibatch] += qpair
+                        lennard_jones[ibatch] += ljpair
                     end
                 end
             end
         end
     end
-    return (332.05382e0 * 4.184) * sum(electrostatic_potential)
+    return (332.05382e0 * 4.184) * sum(electrostatic_potential), sum(lennard_jones)
 end
 
 include("./testing.jl")