Skip to content

Commit

Permalink
a few diff changes, few WIPs
Browse files Browse the repository at this point in the history
  • Loading branch information
samarjeet committed May 30, 2024
1 parent 2433319 commit c355802
Show file tree
Hide file tree
Showing 762 changed files with 183,812 additions and 76 deletions.
13 changes: 9 additions & 4 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.14)
project(chcuda VERSION 0.1.0
LANGUAGES CXX CUDA)

find_package(CUDA REQUIRED)
#find_package(CUDA REQUIRED)

add_subdirectory(external/pybind11)
#find_package(Torch REQUIRED)
Expand Down Expand Up @@ -46,7 +46,7 @@ target_compile_options(chcudalib PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
# --generate-line-info
# -lineinfo
--expt-relaxed-constexpr
-O3
# -O3
>)
target_compile_features(chcudalib PUBLIC cxx_std_20)

Expand All @@ -63,7 +63,7 @@ target_compile_options(chcudadynlib PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
# --generate-line-info
# -lineinfo
--expt-relaxed-constexpr
-O3
# -O3
>)
set_target_properties(chcudadynlib PROPERTIES LINKER_LANGUAGE CUDA)
target_compile_definitions(chcudadynlib PUBLIC USE_TEXTURE_OBJECTS)
Expand All @@ -72,14 +72,19 @@ target_compile_definitions(chcudadynlib PUBLIC USE_TEXTURE_OBJECTS)
find_library( NVTX_LIBRARY
nvToolsExt
PATHS ENV LD_LIBRARY_PATH )
find_library( CUDA_CUFFT_LIBRARIES
cufft
PATHS ENV LD_LIBRARY_PATH )
target_link_libraries(chcudalib ${NVTX_LIBRARY} ${CUDA_CUFFT_LIBRARIES} ${NETCDF_LIBRARIES})
target_link_libraries(chcudalib ${_Python_LIBRARY_RELEASE})
target_link_libraries(chcudadynlib ${NVTX_LIBRARY} ${CUDA_CUFFT_LIBRARIES} ${NETCDF_LIBRARIES})
target_link_libraries(chcudadynlib ${_Python_LIBRARY_RELEASE})
#target_link_libraries(chcudalib external/tinyxml2)
#target_link_libraries(chcudadynlib external/tinyxml2)

include_directories("${CUDA_INCLUDE_DIRS}")
include_directories("${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES}")
#message(STATUS "CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES: ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES}")
#message(STATUS " {NVTX_LIBRARY} {CUDA_CUFFT_LIBRARIES} " ${NVTX_LIBRARY} ${CUDA_CUFFT_LIBRARIES})

#add_subdirectory(fortran_api)

Expand Down
98 changes: 98 additions & 0 deletions examples/soohyung.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
from charmm import apocharmm as ch
import argparse


def p21_bilayer():
prm = ch.CharmmParameters(
[
"../test/data/par_all36_lipid.prm",
"../test/data/toppar_all36_lipid_cholesterol.str",
"../test/data/toppar_water_ions.str",
]
)

filePath = "/v/gscratch/mbs/spark/induced-order/test_p21_apo_charmm/sq2x2/low/"
psf = ch.CharmmPSF(filePath + "init.psf")
# minimized for P21
crd = ch.CharmmCrd(filePath + "step6.6_equilibration.crd")

fm = ch.ForceManager(psf, prm)
fm.setBoxDimensions([144.96, 144.96, 98.65])
fm.setFFTGrid(150, 150, 108)
fm.setKappa(0.34)
fm.setCutoff(14.0) # cutnb of CHARMM
fm.setCtonnb(10.0)
fm.setCtofnb(12.0)
fm.setPeriodicBoundaryCondition(ch.P21) # Only for P21 simulations

ctx = ch.CharmmContext(fm)

ctx.setCoordinates(crd)
ctx.assignVelocitiesAtTemperature(300)

langevinThermostat = ch.LangevinThermostatIntegrator(0.002)
langevinThermostat.setFriction(12.0)
langevinThermostat.setBathTemperature(298.17)
langevinThermostat.setSimulationContext(ctx)

subscriber = ch.DcdSubscriber("out/p21_soohyung_nvt_new.dcd", 1000)
# restartsubscriber = ch.RestartSubscriber("out/p21_soohyung_nvt.res", 10000)

langevinThermostat.subscribe([subscriber]) # , restartsubscriber])

numSteps = int(5e4)
langevinThermostat.propagate(numSteps)

langevinPiston = ch.LangevinPistonIntegrator(0.002)
langevinPiston.setCrystalType(ch.CRYSTAL.TETRAGONAL)
systemMass = psf.getMass()
pistonMass = 0.05 * systemMass

langevinPiston.setPistonMass(
[pistonMass, pistonMass]) # 2 values for x/y and z
langevinPiston.setPistonFriction(20.0)
langevinPiston.setBathTemperature(298.17)
langevinPiston.setSimulationContext(ctx)

"""
restartsubscriber = ch.RestartSubscriber("old.res", 10000) # an already existing restart file
langevinPiston.subscribe([restartsubscriber])
restartsubscriber.readRestart()
"""

baseName = f"out/p21_soohyung_npt_new_{0.002}_{20.0}_{pistonMass}"

dcdSubscriber = ch.DcdSubscriber(f"{baseName}.dcd", 1000)
restartsubscriber = ch.RestartSubscriber(f"{baseName}.res", 10000)
# stateSub = ch.StateSubscriber(f'{baseName}.state', 1000)
# stateSub.setReportFlags("all")

# langevinPiston.subscribe([subscriber, restartsubscriber, stateSub])

langevinPiston.subscribe([dcdSubscriber]) # , restartsubscriber])

nptSteps = int(5e6)
langevinPiston.propagate(nptSteps)


if __name__ == "__main__":
"""
parser = argparse.ArgumentParser(description="P21 bilayer simulation")
parser.add_argument("timeStep", type=float, help="Time step value")
parser.add_argument(
"frictionCoefficient", type=float, help="Friction coefficient value"
)
parser.add_argument(
"pressurePistonMassFactor",
type=float,
help="Pressure piston mass factor : X default value",
)
parser.add_argument(
"temeperaturePistonMassFactor",
type=float,
help="Temperature piston mass factor : X value",
)
args = parser.parse_args()
"""

p21_bilayer()
2 changes: 2 additions & 0 deletions include/CharmmContext.h
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,8 @@ class CharmmContext : public std::enable_shared_from_this<CharmmContext> {
*/
int getDegreesOfFreedom();

// int getMass();

/** @brief Returns the int numDegreesOfFreedom. */
int getNumDegreesOfFreedom();

Expand Down
31 changes: 23 additions & 8 deletions include/CharmmPSF.h
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,12 @@ class Dihedral {
int atom1, atom2, atom3, atom4;
};

class CrossTerm {
public:
int atomi1, atomj1, atomk1, atoml1;
int atomi2, atomj2, atomk2, atoml2;
};

struct InclusionExclusion {
std::vector<int> sizes;
std::vector<int> in14_ex14;
Expand Down Expand Up @@ -74,7 +80,8 @@ class CharmmPSF {
/**
* @brief Increase the mass of the hydrogen atoms to a given amount
* Note that function is not 'repartitioning' hydrogen mass to the neighbor
* atoms, but rather increasing the mass of the hydrogen atoms to a given mass
* atoms, but rather increasing the mass of the hydrogen atoms to a given
* mass
* @param _newHyrogenMass (in a.m.u.)
*/
void setHydrogenMass(double _newHyrogenMass);
Expand All @@ -84,6 +91,7 @@ class CharmmPSF {
int getNumAngles() { return numAngles; }
int getNumDihedrals() { return numDihedrals; }
int getNumImpropers() { return numImpropers; }
int getNumCrossTerms() { return numCrossTerms; }

/** @brief Returns a vector containing all Bond objects */
std::vector<Bond> getBonds() { return bonds; }
Expand All @@ -94,12 +102,14 @@ class CharmmPSF {
/** @brief Returns a vector containing all improper dihedral objects */
std::vector<Dihedral> getImpropers() { return impropers; }

std::vector<CrossTerm> getCrossTerms() { return crossTerms; }

/** @brief Returns a N-sized vector containing all atomic masses */
std::vector<double> getAtomMasses() { return masses; }
/** @brief Returns a N-sized vector containing all atomic charges */
std::vector<double> getAtomCharges() { return charges; }

/** @brief Returns a N-sized string-vector containing all atomic names */
/** @brief Returns a N-sized string-vector containing all atomic names */
std::vector<std::string> getAtomNames() { return atomNames; }
/** @brief Returns a N-sized string-vector containing all atomic types */
std::vector<std::string> getAtomTypes() { return atomTypes; }
Expand All @@ -120,11 +130,14 @@ class CharmmPSF {

/** @brief Computes number of degrees of freedom (3N-6)
*
* Handling constraints is not done here. Should be done by the ForceManager.
* Handling constraints is not done here. Should be done by the
* ForceManager.
*/

int getDegreesOfFreedom();

int getMass();

InclusionExclusion getInclusionExclusionLists();

/** @brief Set charges value */
Expand All @@ -140,15 +153,15 @@ class CharmmPSF {
const std::vector<std::string> &sequence, std::string segmnet);

/**
* @brief Appends the sequence of residues to the PSF structure using the RTF
* object
* @brief Appends the sequence of residues to the PSF structure using the
* RTF object
* @param same as generate
* If @param segment is same already present in the psf, the residues will be
* appeneded to it. Otherwise, a segment in the psf will be created
* If @param segment is same already present in the psf, the residues will
* be appeneded to it. Otherwise, a segment in the psf will be created
*/
void append(const CharmmResidueTopology &rtf,
const std::vector<std::string> &sequence, std::string segmnet);

std::string getOriginalPSFFileName() { return originalPSFFileName; }

private:
Expand All @@ -157,6 +170,7 @@ class CharmmPSF {
int numAngles;
int numDihedrals;
int numImpropers;
int numCrossTerms;
/** @brief N-sized vector containing the atomic masses */
std::vector<double> masses;

Expand All @@ -173,6 +187,7 @@ class CharmmPSF {
std::vector<Angle> angles;
std::vector<Dihedral> dihedrals;
std::vector<Dihedral> impropers;
std::vector<CrossTerm> crossTerms;

CudaContainer<int4> waterMolecules;
CudaContainer<int2> residues;
Expand Down
50 changes: 48 additions & 2 deletions include/CharmmParameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,7 @@ struct BondKey {
return this->atom2 == other.atom2 && this->atom1 == other.atom1;
}
friend std::ostream &operator<<(std::ostream &output, const BondKey &key) {
output << key.atom1 << " " << key.atom2 << " "
<< " ";
output << key.atom1 << " " << key.atom2 << " " << " ";
return output;
}
};
Expand Down Expand Up @@ -136,6 +135,16 @@ struct ImDihedralValues {
}
};

struct CmapKey {
// CmapKey(std::string a1i, std::string a1j, std::string a1k, std::string a1l,
// std::string a2i, std::string a2j, std::string a2k, std::string a2l)
// : atom1i(a1i), atom1j(a1j), atom1k(a1k), atom1l(a1l), atom2i(a2i),
// atom2j(a2j), atom2k(a2k), atom2l(a2l) {}
// std::string atom1i, atom1j, atom1k, atom1l, atom2i, atom2j, atom2k, atom2l;
CmapKey(DihedralKey d1, DihedralKey d2) : dih1(d1), dih2(d2) {}
DihedralKey dih1, dih2;
};

class VdwParameters {
public:
VdwParameters() = default;
Expand All @@ -149,6 +158,39 @@ class VdwParameters {
}
};

class NBFixParameters {
public:
std::string atom1, atom2;
double emin, rmin, emin14, rmin14;
// NBFixParameters(std::string a1, std::string a2, double e, double r,
// double e14, double r14)
// : atom1(a1), atom2(a2), emin(e), rmin(r), emin14(e14), rmin14(r14) {}
// // copy constructor
// NBFixParameters(const NBFixParameters &nbfix) = default;
// NBFixParameters(const NBFixParameters &nbfix)
// : atom1(nbfix.atom1), atom2(nbfix.atom2), emin(nbfix.emin),
// rmin(nbfix.rmin), emin14(nbfix.emin14), rmin14(nbfix.rmin14) {}

// // assignment operator
// NBFixParameters &operator=(const NBFixParameters &nbfix) {
// atom1 = nbfix.atom1;
// atom2 = nbfix.atom2;
// emin = nbfix.emin;
// rmin = nbfix.rmin;
// emin14 = nbfix.emin14;
// rmin14 = nbfix.rmin14;
// return *this;
// }

friend std::ostream &operator<<(std::ostream &output,
const NBFixParameters &nbfix) {
output << "(" << nbfix.atom1 << "," << nbfix.atom2 << "," << nbfix.emin
<< "," << nbfix.rmin << "," << nbfix.emin14 << "," << nbfix.rmin14
<< ")\n";
return output;
}
};

/** @brief Contains bonded interactions parameters and list
* @todo improve doc
*
Expand Down Expand Up @@ -252,6 +294,10 @@ class CharmmParameters {
std::map<DihedralKey, std::vector<DihedralValues>> dihedralParams;
std::map<DihedralKey, ImDihedralValues> improperParams;

std::map<std::tuple<std::string, std::string>, NBFixParameters>
// std::map<BondKey, NBFixParameters>
nbfixParams; // will be filled in the vdw(14)Params

std::map<std::string, VdwParameters> vdwParams;
std::map<std::string, VdwParameters> vdw14Params;
/** @brief original prm file names */
Expand Down
32 changes: 32 additions & 0 deletions include/EDSSubscriber.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
// BEGINLICENSE
//
// This file is part of chcuda, which is distributed under the BSD 3-clause
// license, as described in the LICENSE file in the top level directory of this
// project.
//
// Author: Samarjeet Prasad
//
// ENDLICENSE

#pragma once

#include "ForceManager.h"
#include "Subscriber.h"
#include <fstream>

/** @brief Linked to several ForceManager objects. When reporting, asks each FM
* to compute energy and reports it along with the actual CharmmContext FM. */

/** @brief Reports the enegy of every single ForceManager child */
class EDSSubscriber : public Subscriber {
public:
EDSSubscriber(const std::string &fileName);
EDSSubscriber(const std::string &fileName, int reportFreq);
void update() override;
~EDSSubscriber();

private:
void initialize();
int numFramesWritten;
std::vector<std::shared_ptr<ForceManager>> fmlist;
};
Loading

0 comments on commit c355802

Please sign in to comment.