diff --git a/Jenkinsfile b/Jenkinsfile index 0f40dfc1e2..350b2da092 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -191,56 +191,56 @@ pipeline { } } - stage("Post-test steps") { - parallel { - stage("Publish the manual") { - stages { - stage("Build the manual") { - agent { - docker { - image 'ambermd/lyx:latest' - alwaysPull true - } - } - - steps { - unstash "source" - sh """#!/bin/sh -ex - make docs - cd doc - lyx -batch --export pdf2 CpptrajManual.lyx - lyx -batch --export pdf2 CpptrajDevelopmentGuide.lyx - """ - stash includes: "doc/**", name: "documentation" - } - - post { cleanup { deleteDir() } } - } - stage("Publish the manual") { - agent { label 'linux' } - - steps { - unstash 'documentation' - // Eventually it would be nice to do something better than simply archive - // and make the artifacts available from Jenkins. - archiveArtifacts 'doc/CpptrajManual.pdf,doc/CpptrajDevelopmentGuide.pdf' - // It would be nice to get this in a better place, but for now this - // will suffice. - publishHTML([ - allowMissing: false, - alwaysLinkToLastBuild: false, - keepAll: false, - reportDir: 'doc/html', - reportFiles: 'index.html', - reportName: 'Doxygen Documentation', - reportTitles: '']) - } - - post { cleanup { deleteDir() } } - } - } - } - } - } + //stage("Post-test steps") { + // parallel { + // stage("Publish the manual") { + // stages { + // stage("Build the manual") { + // agent { + // docker { + // image 'ambermd/lyx:latest' + // alwaysPull true + // } + // } + + // steps { + // unstash "source" + // sh """#!/bin/sh -ex + // make docs + // cd doc + // lyx -batch --export pdf2 CpptrajManual.lyx + // lyx -batch --export pdf2 CpptrajDevelopmentGuide.lyx + // """ + // stash includes: "doc/**", name: "documentation" + // } + + // post { cleanup { deleteDir() } } + // } + // stage("Publish the manual") { + // agent { label 'linux' } + + // steps { + // unstash 'documentation' + // // Eventually it would be nice to do something better than simply archive + // // and make the artifacts available from Jenkins. + // archiveArtifacts 'doc/CpptrajManual.pdf,doc/CpptrajDevelopmentGuide.pdf' + // // It would be nice to get this in a better place, but for now this + // // will suffice. + // publishHTML([ + // allowMissing: false, + // alwaysLinkToLastBuild: false, + // keepAll: false, + // reportDir: 'doc/html', + // reportFiles: 'index.html', + // reportName: 'Doxygen Documentation', + // reportTitles: '']) + // } + + // post { cleanup { deleteDir() } } + // } + // } + // } + // } + //} } } diff --git a/README.md b/README.md index 185975b7a8..491fe4fbe2 100644 --- a/README.md +++ b/README.md @@ -53,17 +53,14 @@ specified in the file LICENSE. Documentation ============= The `/doc` subdirectory contains PDF and LyX versions of the CPPTRAJ manual. -The latest version of the manual is compiled each time new code is merged -and is available for download -[here.](https://jenkins.jasonswails.com/job/amber-github/job/cpptraj/job/master/lastSuccessfulBuild/artifact/doc/CpptrajManual.pdf) +The latest version of the manual is available for download +[here.](https://raw.githubusercontent.com/Amber-MD/cpptraj/master/doc/CpptrajManual.pdf) An HTML version can be found [here](https://amber-md.github.io/cpptraj/). There is also limited help for commands in interactive mode via `help []`; `help` with no arguments lists all known commands. - Code documentation generated by doxygen are available -[here.](https://jenkins.jasonswails.com/job/amber-github/job/cpptraj/job/PR-773/Doxygen_20Documentation/) -You can generate the documentation yourself with the command `make docs`. A -limited developers guide is available [here](https://jenkins.jasonswails.com/job/amber-github/job/cpptraj/job/master/lastSuccessfulBuild/artifact/doc/CpptrajDevelopmentGuide.pdf) + Code documentation generated by Doxygen can be generated with the command `make docs`. A +limited developers guide is available [here](https://raw.githubusercontent.com/Amber-MD/cpptraj/master/doc/CpptrajDevelopmentGuide.pdf) and limited HTML-formatted documentation is available [here](https://amber-md.github.io/cpptraj/). diff --git a/configure b/configure index 32b41236b3..3ea8373fb6 100755 --- a/configure +++ b/configure @@ -168,7 +168,7 @@ PNETCDF_OPTS='--disable-fortran --disable-cxx' FFTW_SRCTAR='fftw-3.3.9.tar.gz' FFTW_SRCDIR='fftw-3.3.9' -FFTW_URL="ftp://ftp.fftw.org/pub/fftw/$FFTW_SRCTAR" +FFTW_URL="https://fftw.org/pub/fftw/$FFTW_SRCTAR" FFTW_OPTS='--enable-threads --with-combined-threads --with-pic --disable-fortran' # ----- Variables for downloading MPI ------------ diff --git a/doc/CpptrajDevelopmentGuide.pdf b/doc/CpptrajDevelopmentGuide.pdf new file mode 100644 index 0000000000..2a583ad7c7 Binary files /dev/null and b/doc/CpptrajDevelopmentGuide.pdf differ diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf new file mode 100644 index 0000000000..0ca6f34a70 Binary files /dev/null and b/doc/CpptrajManual.pdf differ diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 2930deb75b..ce198b34a8 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -14247,7 +14247,7 @@ The following topology related commands are available: \begin_layout Standard \align center \begin_inset Tabular - + @@ -14345,6 +14345,26 @@ Print bond info for selected atoms. \begin_inset Text +\begin_layout Plain Layout +bondparminfo +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Print the bond parameter table. +\end_layout + +\end_inset + + + + +\begin_inset Text + \begin_layout Plain Layout change \end_layout @@ -15057,6 +15077,67 @@ If 2 masks are given instead of 1, print info for bonds with first atom in and second atom in . \end_layout +\begin_layout Subsection +bondparminfo +\end_layout + +\begin_layout LyX-Code +bondparminfo [parm | crdset | parmindex <#> | <#>] [out ] +\end_layout + +\begin_deeper +\begin_layout Description +[parm +\begin_inset space ~ +\end_inset + + +\begin_inset space ~ +\end_inset + +| +\begin_inset space ~ +\end_inset + +parmindex +\begin_inset space ~ +\end_inset + +<#> +\begin_inset space ~ +\end_inset + +| +\begin_inset space ~ +\end_inset + +<#>] Name/tag or index of topology. + Default is first loaded topology. +\end_layout + +\begin_layout Description +[out +\begin_inset space ~ +\end_inset + +] File to print to (default STDOUT). +\end_layout + +\end_deeper +\begin_layout Standard +Print the bond parameter table with format: +\end_layout + +\begin_layout LyX-Code +#Idx Rk Req +\end_layout + +\begin_layout Standard +Where Idx is the internal bond parameter index, + Rk is the bond force constant (in kcal/mol*Ang^2), + and Req is the bond equilibrium value (in Ang.). +\end_layout + \begin_layout Subsection change \end_layout @@ -43023,7 +43104,7 @@ vector [] [out [ptrajoutput]] [] [] \end_layout \begin_layout LyX-Code - [magnitude] [ired] [gridset ] + [magnitude] [geom] [ired] [gridset ] [debye] \end_layout \begin_layout LyX-Code @@ -43103,6 +43184,11 @@ ptraj [magnitude] Store the magnitude of the vector with aspect [Mag]. \end_layout +\begin_layout Description +[geom] If specified, + use geometric centers instead of centers of mass. +\end_layout + \begin_layout Description [ired] Mark this vector for subsequent IRED analysis with commands 'matrix ired' and 'ired'. \end_layout @@ -43129,6 +43215,14 @@ ucell[x|y|z] . \end_layout +\begin_layout Description +[debye] ( +\series bold +'dipole' +\series default + vector only) Report dipole vector in units of Debye instead of e-*Ang. +\end_layout + \begin_layout Standard Data Sets Created: \end_layout @@ -43186,8 +43280,13 @@ dipole Store the dipole and center of mass of the atoms specified in \series default . - The vector is not converted to appropriate units, - nor is the value well-defined if the atoms in the mask are not overall charge neutral. + The dipole vector has units of e-*Ang unless +\series bold +'debye' +\series default + is specified for units of Debye. + The center is always stored as simply Ang (since it is just coordinates). + Note that the value may not be well-defined if the atoms in the mask are not overall charge neutral. \end_layout \begin_layout Description diff --git a/src/Action_Angle.cpp b/src/Action_Angle.cpp index 5cd0c61ce8..30f6365045 100644 --- a/src/Action_Angle.cpp +++ b/src/Action_Angle.cpp @@ -62,6 +62,12 @@ Action::RetType Action_Angle::Setup(ActionSetup& setup) { mprintf("Warning: angle: One or more masks contain 0 atoms.\n"); return Action::SKIP; } + // Check if center of mass is possible + if (useMass_) { + if (setup.Top().MaskHasZeroMass(Mask1_)) useMass_ = false; + if (setup.Top().MaskHasZeroMass(Mask2_)) useMass_ = false; + if (setup.Top().MaskHasZeroMass(Mask3_)) useMass_ = false; + } return Action::OK; } diff --git a/src/Action_Dihedral.cpp b/src/Action_Dihedral.cpp index 2ce06fe3e0..db25fde0e6 100644 --- a/src/Action_Dihedral.cpp +++ b/src/Action_Dihedral.cpp @@ -91,7 +91,13 @@ Action::RetType Action_Dihedral::Setup(ActionSetup& setup) { mprintf("Warning: One or more masks have no atoms.\n"); return Action::SKIP; } - + // Check if center of mass is possible + if (useMass_) { + if (setup.Top().MaskHasZeroMass(M1_)) useMass_ = false; + if (setup.Top().MaskHasZeroMass(M2_)) useMass_ = false; + if (setup.Top().MaskHasZeroMass(M3_)) useMass_ = false; + if (setup.Top().MaskHasZeroMass(M4_)) useMass_ = false; + } return Action::OK; } diff --git a/src/Action_Distance.cpp b/src/Action_Distance.cpp index 054497d20a..e915ce5efc 100644 --- a/src/Action_Distance.cpp +++ b/src/Action_Distance.cpp @@ -136,6 +136,13 @@ Action::RetType Action_Distance::Setup(ActionSetup& setup) { else mprintf(", imaging off"); mprintf(".\n"); + // Check if center of mass not possible + if (useMass_) { + if (setup.Top().MaskHasZeroMass( Mask1_ )) useMass_ = false; + if (mode_ == NORMAL) { + if (setup.Top().MaskHasZeroMass( Mask2_ )) useMass_ = false; + } + } return Action::OK; } diff --git a/src/Action_Vector.cpp b/src/Action_Vector.cpp index def85647e7..6d0b3917df 100644 --- a/src/Action_Vector.cpp +++ b/src/Action_Vector.cpp @@ -15,6 +15,8 @@ Action_Vector::Action_Vector() : vcorr_(0), ptrajoutput_(false), needBoxInfo_(false), + useMass_(true), + dipole_in_debye_(false), CurrentParm_(0), outfile_(0) {} @@ -22,7 +24,7 @@ Action_Vector::Action_Vector() : // Action_Vector::Help() void Action_Vector::Help() const { mprintf("\t[] [out [ptrajoutput]] [] []\n" - "\t[magnitude] [ired] [gridset ]\n" + "\t[magnitude] [geom] [ired] [gridset ] [debye]\n" "\t = { mask | minimage | dipole | center | corrplane | \n" "\t box | boxcenter | ucellx | ucelly | ucellz | \n" "\t momentum | principal [x|y|z] | velocity | force }\n" @@ -38,7 +40,10 @@ void Action_Vector::Help() const { " momentum : Store total momentum vector of atoms in (requires velocities).\n" " principal [x|y|z]: X, Y, or Z principal axis vector for atoms in .\n" " velocity : Store velocity of atoms in (requires velocities).\n" - " force : Store force of atoms in (requires forces).\n"); + " force : Store force of atoms in (requires forces).\n" + " If 'magnitude' is specified, also calculate the vector magnitude.\n" + " If 'geom' is specified, use geometric centers, otherwise use center of mass.\n" + " If 'debye' is specified, report 'dipole' vector in Debye instead of e-*Ang.\n"); } // DESTRUCTOR @@ -97,6 +102,7 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d mprinterr("Error: 'ptrajoutput' and 'magnitude' are incompatible.\n"); return Action::ERR; } + dipole_in_debye_ = actionArgs.hasKey("debye"); needBoxInfo_ = false; // Deprecated: corrired, corr, ired if ( actionArgs.hasKey("principal") ) { @@ -209,6 +215,12 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d mprintf("\n"); if (gridSet_ != 0) mprintf("\tExtracting box vectors from grid set '%s'\n", gridSet_->legend()); + if (mode_ == DIPOLE) { + if (dipole_in_debye_) + mprintf("\tDipole vector units will be Debye.\n"); + else + mprintf("\tDipole vector units will be e-*Ang.\n"); + } return Action::OK; } @@ -257,6 +269,11 @@ Action::RetType Action_Vector::Setup(ActionSetup& setup) { return Action::ERR; } } + // Check if center of mass is possible + if (useMass_) { + if (mask_.MaskStringSet() && setup.Top().MaskHasZeroMass( mask_ )) useMass_ = false; + if (mask2_.MaskStringSet() && setup.Top().MaskHasZeroMass( mask2_ )) useMass_ = false; + } CurrentParm_ = setup.TopAddress(); return Action::OK; } @@ -371,10 +388,19 @@ Vec3 Action_Vector::leastSquaresPlane(int n, const double* vcorr) { } // ----------------------------------------------------------------------------- +Vec3 Action_Vector::GetVec(Frame const& currentFrame, AtomMask const& maskIn) +const +{ + if (useMass_) + return currentFrame.VCenterOfMass(maskIn); + else + return currentFrame.VGeometricCenter(maskIn); +} + // Action_Vector::Mask() void Action_Vector::Mask(Frame const& currentFrame) { - Vec3 CXYZ = currentFrame.VCenterOfMass(mask_); - Vec3 VXYZ = currentFrame.VCenterOfMass(mask2_); + Vec3 CXYZ = GetVec(currentFrame, mask_); + Vec3 VXYZ = GetVec(currentFrame, mask2_); VXYZ -= CXYZ; Vec_->AddVxyzo(VXYZ, CXYZ); } @@ -396,6 +422,8 @@ void Action_Vector::Dipole(Frame const& currentFrame) { VXYZ += ( XYZ ); } CXYZ /= total_mass; + if (dipole_in_debye_) + VXYZ /= Constants::DEBYE_EA; Vec_->AddVxyzo( VXYZ, CXYZ ); } @@ -421,7 +449,7 @@ void Action_Vector::Principal(Frame const& currentFrame) { // Action_Vector::CorrPlane() void Action_Vector::CorrPlane(Frame const& currentFrame) { - Vec3 CXYZ = currentFrame.VCenterOfMass(mask_); + Vec3 CXYZ = GetVec(currentFrame, mask_); int idx = 0; for (AtomMask::const_iterator atom = mask_.begin(); atom != mask_.end(); ++atom) @@ -454,8 +482,8 @@ void Action_Vector::BoxLengths(Box const& box) { // Action_Vector::MinImage() void Action_Vector::MinImage(Frame const& frm) { - Vec3 com1 = frm.VCenterOfMass(mask_); - Vec_->AddVxyzo( MinImagedVec(com1, frm.VCenterOfMass(mask2_), frm.BoxCrd().UnitCell(), frm.BoxCrd().FracCell()), com1 ); + Vec3 com1 = GetVec(frm, mask_); + Vec_->AddVxyzo( MinImagedVec(com1, GetVec(frm, mask2_), frm.BoxCrd().UnitCell(), frm.BoxCrd().FracCell()), com1 ); } /// \return The center of selected elements in given array. @@ -476,7 +504,7 @@ static inline Vec3 CalcCenter(const double* xyz, AtomMask const& maskIn) { Action::RetType Action_Vector::DoAction(int frameNum, ActionFrame& frm) { switch ( mode_ ) { case MASK : Mask(frm.Frm()); break; - case CENTER : Vec_->AddVxyz( frm.Frm().VCenterOfMass(mask_) ); break; + case CENTER : Vec_->AddVxyz( GetVec(frm.Frm(), mask_) ); break; case MOMENTUM : Vec_->AddVxyz( frm.Frm().VMomentum(mask_) ); break; case VELOCITY : Vec_->AddVxyz( CalcCenter(frm.Frm().vAddress(), mask_) ); break; case FORCE : Vec_->AddVxyz( CalcCenter(frm.Frm().fAddress(), mask_) ); break; diff --git a/src/Action_Vector.h b/src/Action_Vector.h index 471bc79e72..3af97fdfc1 100644 --- a/src/Action_Vector.h +++ b/src/Action_Vector.h @@ -26,6 +26,8 @@ class Action_Vector : public Action { static double solve_cubic_eq(double,double,double,double); static Vec3 leastSquaresPlane(int,const double*); + /// \return Center of mass or geometric center of atoms in given mask + inline Vec3 GetVec(Frame const&, AtomMask const&) const; void Mask(Frame const&); void Dipole(Frame const&); void Principal(Frame const&); @@ -40,7 +42,9 @@ class Action_Vector : public Action { double* vcorr_; ///< Temp. space for calculating CorrPlane vectorMode mode_; ///< Vector calculation mode bool ptrajoutput_; ///< If true output in ptraj format - bool needBoxInfo_; ///< If true box info required. + bool needBoxInfo_; ///< If true box info required. + bool useMass_; ///< If true, centers are mass-weighted + bool dipole_in_debye_; ///< If true, report dipole vector values in Debye Topology* CurrentParm_; ///< Current topology (for dipole) AtomMask mask_; AtomMask mask2_; diff --git a/src/Atom.cpp b/src/Atom.cpp index 15f8500181..4050c73b5a 100644 --- a/src/Atom.cpp +++ b/src/Atom.cpp @@ -17,7 +17,7 @@ const int Atom::AtomicElementNum_[NUMELEMENTS_] = { 0, 38, 50, 51, 22, 43, 52, 73, 81, 23, 74, 54, 40, 39, 71, 62, - 0 + 0, 0 }; /// Atom names corresponding to AtomicElementType. @@ -36,7 +36,7 @@ const char* Atom::AtomicElementName_[NUMELEMENTS_] = { "??", "SR", "SN", "SB", "TI", "TC", "TE", "TA", "TL", "V", "W", "XE", "ZR", "Y", "LU", "SM", - "XP" + "XP", "DR" }; /** Values taken from 'http://www.webelements.com/' */ @@ -54,7 +54,7 @@ const double Atom::AtomicElementMass_[NUMELEMENTS_] = { 1.0, 87.62, 118.710, 121.760, 47.867, 98, 127.60, 180.94788, 204.3833, 50.9415, 183.84, 131.293, 91.224, 88.90585, 174.9668, 150.36, - 0.0 + 0.0, 0.0 }; /** Values taken from parm10.dat and ion frcmod files as the average Rmin/2 @@ -72,7 +72,7 @@ const double Atom::AtomicElementRadius_[NUMELEMENTS_] = { 1.0, 1.745, 1.000, 1.000, 1.000, 1.000, 1.000, 2.019, 2.419, 1.000, 1.000, 1.000, 1.666, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.850, - 0.000 /* extra point has no radius */ + 0.000, 0.0 /* extra point/Drude has no radius */ }; // CONSTRUCTOR @@ -141,7 +141,9 @@ void Atom::SetMassFromElement() { mass_ = AtomicElementMass_[ element_ ]; } -// CONSTRUCTOR +/** CONSTRUCTOR - Take name, charge, mass, type. Guess element. + * TODO - deprecate this, currently only retained for pytraj. + */ Atom::Atom(NameType const& aname, double charge, double mass, NameType const& atype) : charge_(charge), polar_(0.0), mass_(mass), gb_radius_(0.0), gb_screen_(0.0), aname_(aname), atype_(atype), atype_index_(0), element_(UNKNOWN_ELEMENT), @@ -150,6 +152,18 @@ Atom::Atom(NameType const& aname, double charge, double mass, NameType const& at DetermineElement(0); } +/** CONSTRUCTOR - name, charge, mass, type, element. If element + * is unknown, guess. + */ +Atom::Atom(NameType const& aname, double charge, double mass, NameType const& atype, AtomicElementType elt) : + charge_(charge), polar_(0.0), mass_(mass), gb_radius_(0.0), gb_screen_(0.0), + aname_(aname), atype_(atype), atype_index_(0), element_(elt), + resnum_(0), mol_(0) +{ + if (element_ == UNKNOWN_ELEMENT) + DetermineElement(0); +} + // CONSTRUCTOR Atom::Atom( NameType const& name, double charge, double polar, int atomicnum, double mass, int atidx, NameType const& type, double rad, double screen) : @@ -625,7 +639,13 @@ double Atom::GetBondLength(AtomicElementType atom1, AtomicElementType atom2) { } } else { AtomicElementType e1, e2; - if (atom1 < atom2) { + if (atom1 == EXTRAPT) { + e1 = atom1; + e2 = atom1; + } else if (atom2 == EXTRAPT) { + e1 = atom2; + e2 = atom1; + } else if (atom1 < atom2) { e1 = atom1; e2 = atom2; } else { @@ -700,6 +720,16 @@ double Atom::GetBondLength(AtomicElementType atom1, AtomicElementType atom2) { default: WarnBondLengthDefault(e1,e2,cut); } break; + case EXTRAPT: // Bonds to lone pairs + switch (e2) { + case OXYGEN : cut=0.3; break; + case CARBON : cut=1.4; break; + default: // Different default for lone pairs + cut=0.4; + mprintf("Warning: Bond length not found for %s - %s, using default= %f\n", + AtomicElementName_[e1], AtomicElementName_[e2], cut); + } + break; default: WarnBondLengthDefault(e1,e2,cut); } // END switch(e1) } diff --git a/src/Atom.h b/src/Atom.h index ddc79dad78..03cf6cf666 100644 --- a/src/Atom.h +++ b/src/Atom.h @@ -28,7 +28,7 @@ class Atom { STRONTIUM, TIN, ANTIMONY, TITANIUM, TECHNETIUM, TELLURIUM, // 66 TANTALUM, THALLIUM, VANADIUM, TUNGSTEN, XENON, ZIRCONIUM, // 72 YTTRIUM, LUTETIUM, SAMARIUM, - EXTRAPT + EXTRAPT, DRUDE }; // Constructors and assignment --------------- Atom(); @@ -39,8 +39,10 @@ class Atom { Atom(NameType const&, NameType const&, double); /// Take atom name, type name, and type index. Atom(NameType const&, NameType const&, int); - /// Take atom name, charge, mass, and type name + /// Take atom name, charge, mass, and type name (TODO deprecate) Atom(NameType const&, double, double, NameType const&); + /// Take atom name, charge, mass, type name, and element + Atom(NameType const&, double, double, NameType const&, AtomicElementType); /// Atom name, charge, polarizability, atomic num, mass, type index, type name, gb radius and screen parameters. Atom(NameType const&, double, double, int, double, int, NameType const&, double, double); Atom(const Atom &); @@ -100,7 +102,7 @@ class Atom { /// Set atomic mass from current element void SetMassFromElement(); protected: - static const size_t NUMELEMENTS_ = 77; + static const size_t NUMELEMENTS_ = 78; private: static CPPTRAJ_EXPORT const int AtomicElementNum_[]; static CPPTRAJ_EXPORT const char* AtomicElementName_[]; diff --git a/src/Command.cpp b/src/Command.cpp index d9cf5ecd89..545514dea5 100644 --- a/src/Command.cpp +++ b/src/Command.cpp @@ -297,6 +297,7 @@ void Command::Init() { Command::AddCmd( new Exec_AngleInfo(), Cmd::EXE, 3, "angles", "angleinfo", "printangles" ); Command::AddCmd( new Exec_AtomInfo(), Cmd::EXE, 3, "atoms", "atominfo", "printatoms" ); Command::AddCmd( new Exec_BondInfo(), Cmd::EXE, 3, "bonds", "bondinfo", "printbonds" ); + Command::AddCmd( new Exec_BondParmInfo(), Cmd::EXE, 1, "bondparminfo" ); Command::AddCmd( new Exec_Change(), Cmd::EXE, 1, "change" ); Command::AddCmd( new Exec_ChargeInfo(), Cmd::EXE, 1, "charge" ); Command::AddCmd( new Exec_CompareTop(), Cmd::EXE, 1, "comparetop" ); diff --git a/src/Exec_Top.cpp b/src/Exec_Top.cpp index 253e7ede25..01c22beba7 100644 --- a/src/Exec_Top.cpp +++ b/src/Exec_Top.cpp @@ -82,6 +82,19 @@ Exec::RetType Exec_UBInfo::Execute(CpptrajState& State, ArgList& argIn) { return CpptrajState::OK; } +// ----------------------------------------------------------------------------- +void Exec_BondParmInfo::Help() const { + mprintf("\t[%s] [out ]\n", DataSetList::TopIdxArgs); + mprintf(" Print bond parameter table.\n"); +} + +Exec::RetType Exec_BondParmInfo::Execute(CpptrajState& State, ArgList& argIn) { + TopInfo info; + if (CommonSetup(info, State, argIn, "Bond info")) return CpptrajState::ERR; + if (info.PrintBondParm()) return CpptrajState::ERR; + return CpptrajState::OK; +} + // ----------------------------------------------------------------------------- void Exec_AngleInfo::Help() const { mprintf("\t[%s] [] [ ]\n\t[out ]\n", DataSetList::TopIdxArgs); diff --git a/src/Exec_Top.h b/src/Exec_Top.h index cfcac9a95a..f79ef9ca06 100644 --- a/src/Exec_Top.h +++ b/src/Exec_Top.h @@ -38,6 +38,14 @@ class Exec_UBInfo : public Exec { DispatchObject* Alloc() const { return (DispatchObject*)new Exec_UBInfo(); } RetType Execute(CpptrajState&, ArgList&); }; +/// Print bond parameters +class Exec_BondParmInfo : public Exec { + public: + Exec_BondParmInfo() : Exec(PARM) {} + void Help() const; + DispatchObject* Alloc() const { return (DispatchObject*)new Exec_BondParmInfo(); } + RetType Execute(CpptrajState&, ArgList&); +}; /// Print angle info for atoms in mask. class Exec_AngleInfo : public Exec { public: diff --git a/src/Parm_CharmmPsf.cpp b/src/Parm_CharmmPsf.cpp index 6464acfac1..1abbc424cd 100644 --- a/src/Parm_CharmmPsf.cpp +++ b/src/Parm_CharmmPsf.cpp @@ -129,8 +129,153 @@ int Parm_CharmmPsf::ReadDihedrals(BufferedLine& infile, int ndihedral, const cha return 0; } +// ----- Charmm Lone Pair ------------------------ +class Parm_CharmmPsf::LonePair { + public: + //LonePair() : numSupport_(0), idx_(-1), dist_(0), ang_(0), dih_(0) {} + LonePair(int n, int i, const char* t, double dist, double ang, double dih) : + nat_(n), idx_(i), type_(t), dist_(dist), ang_(ang), dih_(dih) {} + + int Nat() const { return nat_; } + int Idx() const { return idx_; } + const char* lptype() const { return type_.c_str(); } + double Dist() const { return dist_; } + double Ang() const { return ang_; } + double Dih() const { return dih_; } + + void Print() const { mprintf(" %i %i %s %f %f %f", Nat(), Idx()+1, lptype(), Dist(), Ang(), Dih()); } + private: + int nat_; ///< Number of support atoms + int idx_; ///< Index of support atoms in support (host) array + std::string type_; ///< Lone pair type TODO make char? + double dist_; ///< Distance in angstroms + double ang_; ///< Angle in degrees + double dih_; ///< Dihedral in degrees +}; +// ----------------------------------------------- + +/** Read lone pairs. */ +int Parm_CharmmPsf::ReadLonePairs(BufferedLine& infile, int numlp, int numlph, Topology& parmOut) +const +{ + std::vector LPs; + LPs.reserve( numlp ); + // Read the lone pairs + const char* buffer = 0; + char type[9]; + type[8] = '\0'; + int num = 0; + int idx = -1; + double dist = 0; + double ang = 0; + double dih = 0; + for (int ilp = 0; ilp != numlp; ilp++) { + if ( (buffer = infile.Line()) == 0) { + mprinterr("Error: Reading lone pair %i on line %i\n", ilp+1, infile.LineNumber()); + return 1; + } + //mprintf("DEBUG: LP: %s\n", buffer); + num = 0; + idx = -1; + dist = 0; + ang = 0; + dih = 0; + type[0] = '\0'; + int nscan = sscanf(buffer, "%i %i %s %lf %lf %lf", &num, &idx, type, &dist, &ang, &dih); + if (nscan < 4) { + mprinterr("Error: Only read %i values for lone pair line, expected at least 4.\n", nscan); + mprinterr("Error: Line %i: %s\n", infile.LineNumber(), buffer); + return 1; + } + LPs.push_back( LonePair(num, idx-1, type, dist, ang, dih) ); + } + // Read the lone pair atoms + int lpat[8]; + int nlines = numlph / 8; + if ( (numlph%8) > 0 ) + nlines++; + std::vector LPatoms; + LPatoms.reserve( numlph ); + for (int ilp = 0; ilp != nlines; ilp++) { + if ( (buffer = infile.Line()) == 0) { + mprinterr("Error: Reading lone pair atoms on line %i\n", infile.LineNumber()); + return 1; + } + int nscan = sscanf(buffer, "%i %i %i %i %i %i %i %i", + lpat, lpat+1, lpat+2, lpat+3, lpat+4, lpat+5, lpat+6, lpat+7); + for (int ii = 0; ii != nscan; ii++) + LPatoms.push_back( lpat[ii]-1 ); + } + if ( (unsigned int)numlph != LPatoms.size() ) { + mprinterr("Error: Number of lone pair atoms read (%zu) != expected (%i)\n", LPatoms.size(), numlph); + return 1; + } + // Add bonds for lone pairs TODO angles and dihedrals? + for (std::vector::const_iterator it = LPs.begin(); it != LPs.end(); ++it) { + // Sanity check + int lpAtomIdx = LPatoms[it->Idx()]; + int bondedAtomIdx = LPatoms[it->Idx()+1]; + if ( parmOut[lpAtomIdx].Element() != Atom::EXTRAPT ) { + mprintf("Warning: PSF defines atom %s as lone pair, but it does not appear to be one.\n", + parmOut.AtomMaskName( lpAtomIdx ).c_str()); + } + // Bond + if (params_.BP().empty()) { + // Just add the bond to lone pair + parmOut.AddBond( lpAtomIdx, bondedAtomIdx ); + } else { + // Add the distance as a parameter FIXME is Rk=0 ok? + // FIXME what does a negative distance mean? + double req = it->Dist(); + if (req < 0) req = -req; + parmOut.AddBond( lpAtomIdx, bondedAtomIdx, BondParmType(0, req) ); + } + // Angle + if (it->Nat() > 1) { + int angleAtomIdx = LPatoms[it->Idx()+2]; + if (params_.AP().empty()) { + parmOut.AddAngle( lpAtomIdx, bondedAtomIdx, angleAtomIdx ); + } else { + // Add the angle as a parameter (in radians) FIXME is Tk=0 ok? + double teq = it->Ang() * Constants::DEGRAD; + parmOut.AddAngle( lpAtomIdx, bondedAtomIdx, angleAtomIdx, AngleParmType(0, teq) ); + } + } + // Dihedral + if (it->Nat() > 2) { + int angleAtomIdx = LPatoms[it->Idx()+2]; + int dihAtomIdx = LPatoms[it->Idx()+3]; + if (params_.DP().empty()) { + parmOut.AddDihedral( lpAtomIdx, bondedAtomIdx, angleAtomIdx, dihAtomIdx ); + } else { + // Add the dihedral as a parameter (use radians) FIXME Pk, Pn, etc? + DihedralParmType dpt; + dpt.SetPhase( it->Dih() * Constants::DEGRAD ); + // FIXME OK for normal type for lone pair? + parmOut.AddDihedral( DihedralType(lpAtomIdx, bondedAtomIdx, angleAtomIdx, dihAtomIdx, DihedralType::NORMAL), dpt ); + } + } + } + // DEBUG - Print out lone pair information + if (debug_ > 1) { + for (std::vector::const_iterator it = LPs.begin(); it != LPs.end(); ++it) { + mprintf("DEBUG: LP:"); + it->Print(); + mprintf(" Atoms:"); + int jj = it->Idx(); + // Nat() is the number of support atoms; +1 for the lone pair atom itself + for (int ii = 0; ii <= it->Nat(); ii++, jj++) + mprintf(" %i", LPatoms[jj]+1); + mprintf("\n"); + } + } + + return 0; +} + const unsigned int Parm_CharmmPsf::ChmStrMax_ = 9; +/** Extract residue number and alternatively the insertion code. */ int Parm_CharmmPsf::ParseResID(char& psficode, const char* psfresid) { char buf[ChmStrMax_]; @@ -163,7 +308,7 @@ int Parm_CharmmPsf::ParseResID(char& psficode, const char* psfresid) * Mask selection requires natom, nres, names, resnames, resnums. */ int Parm_CharmmPsf::ReadParm(FileName const& fname, Topology &parmOut) { - const size_t TAGSIZE = 10; + const size_t TAGSIZE = 11; char tag[TAGSIZE]; tag[0]='\0'; @@ -251,6 +396,12 @@ int Parm_CharmmPsf::ReadParm(FileName const& fname, Topology &parmOut) { mprintf("Warning: During read of PSF atoms, expected to read 6 columns, got %i (line %i)\n", nread, infile.LineNumber()); } } + // Determine if this is a Drude particle + Atom::AtomicElementType psfElt = Atom::UNKNOWN_ELEMENT; + if (psftype[0] == 'D') { + //mprintf("DEBUG: Potential Drude particle: %s %s\n", psfname, psftype); + psfElt = Atom::DRUDE; + } // Extract residue number and alternatively insertion code. int psfresnum = ParseResID(psficode, psfresid); //mprintf("DEBUG: resnum %10i icode %c\n", psfresnum, psficode); @@ -269,7 +420,7 @@ int Parm_CharmmPsf::ReadParm(FileName const& fname, Topology &parmOut) { } }*/ atomTypes.AddParm( TypeNameHolder(NameType(psftype)), AtomType(psfmass), false ); - Atom chmAtom( psfname, psfcharge, psfmass, psftype ); + Atom chmAtom( psfname, psfcharge, psfmass, psftype, psfElt ); parmOut.AddTopAtom( chmAtom, Residue(psfresname, psfresnum, ' ', std::string(segmentID)) ); } // END loop over atoms // Advance to !NBOND @@ -366,6 +517,20 @@ int Parm_CharmmPsf::ReadParm(FileName const& fname, Topology &parmOut) { if (ReadDihedrals(infile, nimproper, "improper", parmOut)) return 1; } else mprintf("Warning: PSF has no impropers.\n"); + // SKIPPING NDON, NACC, NNB, NGRP + // Advance to <# lone pairs> <# lone pair hosts> NUMLP NUMLPH + int numlp = -1; + int numlph = -1; + while (strncmp(tag, "!NUMLP", 6) !=0) { + const char* buffer = infile.Line(); + if ( buffer == 0 ) break; + sscanf(buffer,"%i %i %10s", &numlp, &numlph, tag); + } + if (numlp > -1) { + if (debug_ > 0) mprintf("DEBUG: PSF contains %i lone pairs, %i lone pair hosts.\n", numlp, numlph); + if (ReadLonePairs(infile, numlp, numlph, parmOut)) return 1; + } + mprintf("\tPSF contains %i atoms, %i residues.\n", parmOut.Natom(), parmOut.Nres()); infile.CloseFile(); diff --git a/src/Parm_CharmmPsf.h b/src/Parm_CharmmPsf.h index ea37fcc1c9..0b5ed474e9 100644 --- a/src/Parm_CharmmPsf.h +++ b/src/Parm_CharmmPsf.h @@ -15,11 +15,14 @@ class Parm_CharmmPsf : public ParmIO { static void WriteHelp(); int processWriteArgs(ArgList&); private: + class LonePair; + static const unsigned int ChmStrMax_; //static inline int FindTag(char*, const char*, int, BufferedLine&); static inline int ParseResID(char&, const char*); static inline int FindTag(char*, const char*, BufferedLine&); int ReadDihedrals(BufferedLine&, int, const char*, Topology&) const; + int ReadLonePairs(BufferedLine&, int, int, Topology&) const; inline void WriteSectionHeader(CpptrajFile&, const char*, int) const; diff --git a/src/Parm_Gromacs.cpp b/src/Parm_Gromacs.cpp index 18f6551a73..3ff8981f77 100644 --- a/src/Parm_Gromacs.cpp +++ b/src/Parm_Gromacs.cpp @@ -414,7 +414,7 @@ int Parm_Gromacs::ReadParm(FileName const& fname, Topology &TopIn) { for (AtomArray::const_iterator atom = Mol.begin(); atom != Mol.end(); ++atom) { if (atom->mass_ > -1.0) - TopIn.AddTopAtom( Atom( atom->aname_, atom->charge_, atom->mass_, atom->atype_), + TopIn.AddTopAtom( Atom( atom->aname_, atom->charge_, atom->mass_, atom->atype_, Atom::UNKNOWN_ELEMENT), Residue(atom->rname_, atom->rnum_ + resoffset, ' ', "") ); else TopIn.AddTopAtom( Atom( atom->aname_, atom->atype_, atom->charge_ ), diff --git a/src/TopInfo.cpp b/src/TopInfo.cpp index 472d7cf8de..123b0a4b27 100644 --- a/src/TopInfo.cpp +++ b/src/TopInfo.cpp @@ -540,6 +540,14 @@ int TopInfo::PrintBondInfo(std::string const& mask1exp, std::string const& mask2 return 0; } +/** Print bond parameters */ +int TopInfo::PrintBondParm() const { + outfile_->Printf("%-10s %8s %8s\n", "#Idx", "Rk", "Req"); + for (BondParmArray::const_iterator bp = parm_->BondParm().begin(); bp != parm_->BondParm().end(); ++bp) + outfile_->Printf("%-10li %8.3f %8.3f\n", bp - parm_->BondParm().begin(), bp->Rk(), bp->Req()); + return 0; +} + // TopInfo::PrintAngles() void TopInfo::PrintAngles(AngleArray const& aarray, AngleParmArray const& angleparm, CharMask const& mask1, CharMask const& mask2, CharMask const& mask3, diff --git a/src/TopInfo.h b/src/TopInfo.h index fe17088ac8..9ac6119278 100644 --- a/src/TopInfo.h +++ b/src/TopInfo.h @@ -24,6 +24,7 @@ class TopInfo { int PrintMoleculeInfo(std::string const&) const; int PrintShortMolInfo(std::string const&) const; int PrintBondInfo(std::string const&, std::string const&, bool) const; + int PrintBondParm() const; int PrintAngleInfo(std::string const&, std::string const&, std::string const&) const; int PrintDihedralInfo(std::string const&, std::string const&, std::string const&, std::string const&,bool) const; diff --git a/src/Topology.cpp b/src/Topology.cpp index 6a5599b6bd..678cc0fd01 100644 --- a/src/Topology.cpp +++ b/src/Topology.cpp @@ -759,6 +759,7 @@ void Topology::AddBond(int atom1, int atom2, BondParmType const& BPin) { pidx = (int)bondparm_.size(); bondparm_.push_back( BPin ); } + //mprintf("DEBUG: Add bond %i - %i Rk=%f Req=%f (%i)\n", atom1+1, atom2+1, BPin.Rk(), BPin.Req(), pidx); AddBond( atom1, atom2, pidx ); } @@ -825,17 +826,43 @@ int Topology::RemoveBond(int atom1, int atom2) * For bonds to H always insert the H second. */ void Topology::AddBond(int atom1, int atom2, int pidxIn) { + //mprintf("DEBUG: Enter AddBond(%i, %i, %i)\n", atom1+1, atom2+1, pidxIn); // Check if atoms are out of range. if (WarnOutOfRange(atoms_.size(), atom1, "bond")) return; if (WarnOutOfRange(atoms_.size(), atom2, "bond")) return; + bool a1H = (atoms_[atom1].Element() == Atom::HYDROGEN); + bool a2H = (atoms_[atom2].Element() == Atom::HYDROGEN); // Check for duplicate bond for (Atom::bond_iterator ba = atoms_[atom1].bondbegin(); ba != atoms_[atom1].bondend(); ++ba) + { if ( *ba == atom2 ) { if (debug_ > 0) mprintf("Warning: Bond between atoms %i and %i already exists.\n", atom1+1, atom2+1); + // If the bond already exists but does not yet have a parameter, update + // the parameter index before exiting. + if (pidxIn > -1 && pidxIn < (int)bondparm_.size()) { + BondArray* BONDS = 0; + if (a1H || a2H) + BONDS = &bondsh_; + else + BONDS= &bonds_; + for (BondArray::iterator it = BONDS->begin(); it != BONDS->end(); ++it) { + if ( (it->A1() == atom1 && it->A2() == atom2) || + (it->A1() == atom2 && it->A2() == atom1) ) + { + mprintf("DEBUG: Existing bond found. Existing Idx %i\n", it->Idx()); + if (it->Idx() < 0) { + mprintf("DEBUG: Adding bond parameter index %i for existing bond.\n", pidxIn); + it->SetIdx( pidxIn ); + } + break; + } + } // END loop over target bond array + } return; } + } // END check for duplicate bond. // Check if parm index is out of range; int pidx; if (pidxIn < (int)bondparm_.size()) @@ -844,9 +871,7 @@ void Topology::AddBond(int atom1, int atom2, int pidxIn) { mprintf("Warning: No bond parameters for index %i\n", pidxIn); pidx = -1; } - bool a1H = (atoms_[atom1].Element() == Atom::HYDROGEN); - bool a2H = (atoms_[atom2].Element() == Atom::HYDROGEN); - //mprintf("\t\t\tAdding bond %i to %i (isH=%i)\n",atom1+1,atom2+1,(int)isH); + //mprintf("\t\t\tAdding bond %i to %i (a1H=%i a2H=%i) idx=%i\n",atom1+1,atom2+1,(int)a1H,(int)a2H,pidx); // Update bonds arrays if (a1H || a2H) { if (a1H) @@ -1397,6 +1422,19 @@ std::vector Topology::MolnumsSelectedBy(AtomMask const& mask) const { return tmp; } +/** \return True if total mass of selected atoms is zero. */ +bool Topology::MaskHasZeroMass(AtomMask const& mask) const { + double totalMass = 0; + for (AtomMask::const_iterator it = mask.begin(); it != mask.end(); ++it) + totalMass += atoms_[*it].Mass(); + if (totalMass == 0) { + mprintf("Warning: The total mass of atoms in mask '%s' is zero; cannot use mass-weighting.\n", + mask.MaskString()); + return true; + } + return false; +} + // ----------------------------------------------------------------------------- int Topology::scale_dihedral_K(DihedralArray& dihedrals, CharMask const& Mask, double scale_factor, bool useAll) diff --git a/src/Topology.h b/src/Topology.h index 8dffa3a128..b13a9d0e92 100644 --- a/src/Topology.h +++ b/src/Topology.h @@ -229,6 +229,8 @@ class Topology { std::vector ResnumsSelectedBy(AtomMask const&) const; /// \return Array of molecule numbers selected by given atom mask std::vector MolnumsSelectedBy(AtomMask const&) const; + /// \return True if the total mass of selected atoms is zero + bool MaskHasZeroMass(AtomMask const&) const; // ----- Topology modification routines ------ int ScaleDihedralK(double, std::string const&, bool); /// Strip atoms outside given mask, do not keep parameters. diff --git a/src/Version.h b/src/Version.h index d1e831b38d..a63722cdb3 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V6.28.0" +#define CPPTRAJ_INTERNAL_VERSION "V6.29.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/test/Test_Charmm/RunTest.sh b/test/Test_Charmm/RunTest.sh index 529f311f1a..6cf385a2de 100755 --- a/test/Test_Charmm/RunTest.sh +++ b/test/Test_Charmm/RunTest.sh @@ -4,7 +4,7 @@ CleanFiles charmm.in test.ala3.pdb.? test.ala3.pdb.10 first.ala3.crd \ test.psf test.ala3.dcd second.ala3.crd strip.chamber.parm7 \ - run0.res_0.mol2 cpptraj.psf cpptraj.cor + run0.res_0.mol2 cpptraj.psf cpptraj.cor ala3.drude.mol2 TESTNAME='Charmm DCD tests' Requires maxthreads 10 @@ -88,6 +88,18 @@ EOF DoTest cpptraj.cor.save cpptraj.cor fi +UNITNAME='CHARMM Drude PSF test' +CheckFor maxthreads 1 +if [ $? -eq 0 ] ; then + cat > charmm.in <MOLECULE +Cpptraj Generated mol2 file. + 70 69 3 0 0 +SMALL +USER_CHARGES + + +@ATOM + 1 CAY -6.3044 -1.9238 0.7786 CD33C 1 ALA 2.286200 + 2 DCAY -6.3039 -1.9240 0.7783 DRUD 1 ALA -2.388200 + 3 HY1 -6.6619 -2.4750 -0.1173 HDA3A 1 ALA 0.048000 + 4 HY2 -6.5481 -0.8465 0.6584 HDA3A 1 ALA 0.048000 + 5 HY3 -6.8358 -2.3133 1.6731 HDA3A 1 ALA 0.048000 + 6 CY -4.8347 -2.0958 0.9489 CD2O1A 1 ALA 1.922700 + 7 DCY -4.8347 -2.0958 0.9489 DRUD 1 ALA -1.425700 + 8 OY -4.3263 -2.0833 2.0675 OD2C1A 1 ALA 1.400200 + 9 DOY -4.3263 -2.0833 2.0675 DRUD 1 ALA -1.400200 + 10 LPY1 -4.5946 -2.0411 2.1948 LPDO1 1 ALA -0.312000 + 11 LPY2 -4.0535 -2.1254 1.9498 LPDO1 1 ALA -0.227000 + 12 N -4.1051 -2.2579 -0.1619 ND2A2 1 ALA 1.959400 + 13 DN -4.1127 -2.2578 -0.1774 DRUD 1 ALA -2.365400 + 14 HN -4.5611 -2.2541 -1.0854 HDP1A 1 ALA 0.296000 + 15 CA -2.6787 -2.3824 -0.2023 CD31C 1 ALA 1.961700 + 16 DCA -2.6744 -2.3877 -0.1902 DRUD 1 ALA -1.791700 + 17 HA -2.3370 -2.8079 0.7655 HDA1A 1 ALA -0.017000 + 18 C -2.0536 -1.0209 -0.3631 CD2O1A 1 ALA 1.860300 + 19 DC -2.0536 -1.0209 -0.3631 DRUD 1 ALA -1.365300 + 20 O -2.6721 -0.1051 -0.9057 OD2C1A 1 ALA 1.320500 + 21 DO -2.6721 -0.1051 -0.9057 DRUD 1 ALA -1.320500 + 22 LPOA -2.9090 -0.2839 -0.9493 LPDO1 1 ALA -0.311000 + 23 LPOB -2.4404 0.0815 -0.8667 LPDO1 1 ALA -0.227000 + 24 CB -2.1995 -3.3012 -1.3384 CD33A 1 ALA 2.231800 + 25 DCB -2.1993 -3.3014 -1.3381 DRUD 1 ALA -2.330800 + 26 HB1 -2.6503 -4.3104 -1.2260 HDA3A 1 ALA 0.033000 + 27 HB2 -2.5066 -2.8917 -2.3244 HDA3A 1 ALA 0.033000 + 28 HB3 -1.0944 -3.4150 -1.3275 HDA3A 1 ALA 0.033000 + 29 N -0.8076 -0.8432 0.0981 ND2A2 2 ALA 1.959400 + 30 DN -0.7993 -0.8552 0.1050 DRUD 2 ALA -2.365400 + 31 HN -0.2754 -1.6073 0.5382 HDP1A 2 ALA 0.296000 + 32 CA -0.1051 0.3983 0.0418 CD31C 2 ALA 1.961700 + 33 DCA -0.1077 0.4043 0.0289 DRUD 2 ALA -1.791700 + 34 HA -0.3016 0.8620 -0.9485 HDA1A 2 ALA -0.017000 + 35 C 1.3638 0.1099 0.1333 CD2O1A 2 ALA 1.860300 + 36 DC 1.3638 0.1099 0.1333 DRUD 2 ALA -1.365300 + 37 O 1.7609 -0.9587 0.5945 OD2C1A 2 ALA 1.320500 + 38 DO 1.7609 -0.9587 0.5945 DRUD 2 ALA -1.320500 + 39 LPOA 1.4788 -1.0494 0.6413 LPDO1 2 ALA -0.311000 + 40 LPOB 2.0464 -0.8771 0.5516 LPDO1 2 ALA -0.227000 + 41 CB -0.5000 1.3629 1.1742 CD33A 2 ALA 2.231800 + 42 DCB -0.4998 1.3634 1.1740 DRUD 2 ALA -2.330800 + 43 HB1 -1.5869 1.5862 1.1196 HDA3A 2 ALA 0.033000 + 44 HB2 -0.2933 0.8997 2.1627 HDA3A 2 ALA 0.033000 + 45 HB3 0.0538 2.3236 1.1045 HDA3A 2 ALA 0.033000 + 46 N 2.2109 1.0553 -0.3021 ND2A2 3 ALA 1.959400 + 47 DN 2.2046 1.0701 -0.3097 DRUD 3 ALA -2.365400 + 48 HN 1.8477 1.9129 -0.7420 HDP1A 3 ALA 0.296000 + 49 CA 3.6377 0.9275 -0.2916 CD31C 3 ALA 1.961700 + 50 DCA 3.6415 0.9208 -0.2794 DRUD 3 ALA -1.791700 + 51 HA 3.9308 0.4061 0.6446 HDA1A 3 ALA -0.017000 + 52 C 4.2992 2.2847 -0.2692 CD2O1A 3 ALA 1.860300 + 53 DC 4.2992 2.2847 -0.2692 DRUD 3 ALA -1.365300 + 54 O 5.5006 2.4051 -0.5048 OD2C1A 3 ALA 1.320500 + 55 DO 5.5006 2.4051 -0.5048 DRUD 3 ALA -1.320500 + 56 LPOA 5.5277 2.1091 -0.5454 LPDO1 3 ALA -0.311000 + 57 LPOB 5.4837 2.7021 -0.4663 LPDO1 3 ALA -0.227000 + 58 CB 4.1513 0.1355 -1.5068 CD33A 3 ALA 2.231800 + 59 DCB 4.1512 0.1354 -1.5071 DRUD 3 ALA -2.330800 + 60 HB1 3.6956 -0.8776 -1.5217 HDA3A 3 ALA 0.033000 + 61 HB2 3.8755 0.6536 -2.4501 HDA3A 3 ALA 0.033000 + 62 HB3 5.2552 0.0151 -1.4740 HDA3A 3 ALA 0.033000 + 63 NT 3.5443 3.3581 0.0100 ND2A2 3 ALA 2.036300 + 64 DNT 3.5274 3.3565 0.0144 DRUD 3 ALA -2.418300 + 65 HNT 2.5532 3.2607 0.2727 HDP1A 3 ALA 0.272000 + 66 CAT 4.0636 4.6931 0.0034 CD33C 3 ALA 2.166600 + 67 DCAT 4.0627 4.6921 0.0039 DRUD 3 ALA -2.221600 + 68 HTC1 4.6262 4.8865 -0.9349 HDA3A 3 ALA 0.055000 + 69 HTC2 3.2382 5.4332 0.0754 HDA3A 3 ALA 0.055000 + 70 HTC3 4.7479 4.8390 0.8664 HDA3A 3 ALA 0.055000 +@BOND + 1 1 2 1 + 2 6 7 1 + 3 8 9 1 + 4 12 6 1 + 5 1 6 1 + 6 6 8 1 + 7 8 10 1 + 8 8 11 1 + 9 12 13 1 + 10 15 16 1 + 11 18 19 1 + 12 20 21 1 + 13 24 25 1 + 14 12 15 1 + 15 15 18 1 + 16 18 29 1 + 17 18 20 1 + 18 15 24 1 + 19 20 22 1 + 20 20 23 1 + 21 29 30 1 + 22 32 33 1 + 23 35 36 1 + 24 37 38 1 + 25 41 42 1 + 26 29 32 1 + 27 32 35 1 + 28 35 46 1 + 29 35 37 1 + 30 32 41 1 + 31 37 39 1 + 32 37 40 1 + 33 46 47 1 + 34 49 50 1 + 35 52 53 1 + 36 54 55 1 + 37 58 59 1 + 38 46 49 1 + 39 49 52 1 + 40 52 54 1 + 41 49 58 1 + 42 54 56 1 + 43 54 57 1 + 44 63 64 1 + 45 66 67 1 + 46 52 63 1 + 47 63 66 1 + 48 1 3 1 + 49 1 4 1 + 50 1 5 1 + 51 15 17 1 + 52 12 14 1 + 53 24 26 1 + 54 24 27 1 + 55 24 28 1 + 56 32 34 1 + 57 29 31 1 + 58 41 43 1 + 59 41 44 1 + 60 41 45 1 + 61 49 51 1 + 62 46 48 1 + 63 58 60 1 + 64 58 61 1 + 65 58 62 1 + 66 63 65 1 + 67 66 68 1 + 68 66 69 1 + 69 66 70 1 +@SUBSTRUCTURE + 1 ALA 1 **** 0 **** **** + 2 ALA 29 **** 0 **** **** + 3 ALA 46 **** 0 **** **** diff --git a/test/Test_Charmm/ala3.drude.psf b/test/Test_Charmm/ala3.drude.psf new file mode 100644 index 0000000000..4f8206e3a0 --- /dev/null +++ b/test/Test_Charmm/ala3.drude.psf @@ -0,0 +1,249 @@ +PSF EXT CMAP DRUDE XPLOR AUTOG + + 4 !NTITLE +* GENERATED BY CHARMM-GUI (HTTP://WWW.CHARMM-GUI.ORG) V3.7 ON AUG, 29. 2024. JOB +* READ CRD FILES AND GENERATE DRUDE PSF +* NB: WATER MUST BE LAST SEGMENT GENERATED DO TO NOANGLE NODIHEDRAL +* DATE: 8/29/24 9:12:20 CREATED BY USER: apache + + 70 !NATOM + 1 ALA3 1 ALA CAY CD33C 2.28620 11.6110 0 -1.89400 2.12200 + 2 ALA3 1 ALA DCAY DRUD -2.38820 0.400000 0 0.00000 0.00000 + 3 ALA3 1 ALA HY1 HDA3A 0.480000E-01 1.00800 0 0.00000 0.00000 + 4 ALA3 1 ALA HY2 HDA3A 0.480000E-01 1.00800 0 0.00000 0.00000 + 5 ALA3 1 ALA HY3 HDA3A 0.480000E-01 1.00800 0 0.00000 0.00000 + 6 ALA3 1 ALA CY CD2O1A 1.92270 11.6110 0 -0.675000 0.295000 + 7 ALA3 1 ALA DCY DRUD -1.42570 0.400000 0 0.00000 0.00000 + 8 ALA3 1 ALA OY OD2C1A 1.40020 15.5994 0 -0.651000 0.310000 + 9 ALA3 1 ALA DOY DRUD -1.40020 0.400000 0 0.00000 0.00000 + 10 ALA3 1 ALA LPY1 LPDO1 -0.312000 0.00000 -1 0.00000 0.00000 + 11 ALA3 1 ALA LPY2 LPDO1 -0.227000 0.00000 -1 0.00000 0.00000 + 12 ALA3 1 ALA N ND2A2 1.95940 13.6070 0 -1.85800 0.126000 + 13 ALA3 1 ALA DN DRUD -2.36540 0.400000 0 0.00000 0.00000 + 14 ALA3 1 ALA HN HDP1A 0.296000 1.00800 0 0.00000 0.00000 + 15 ALA3 1 ALA CA CD31C 1.96170 11.6110 0 -1.06600 1.14100 + 16 ALA3 1 ALA DCA DRUD -1.79170 0.400000 0 0.00000 0.00000 + 17 ALA3 1 ALA HA HDA1A -0.170000E-01 1.00800 0 0.00000 0.00000 + 18 ALA3 1 ALA C CD2O1A 1.86030 11.6110 0 -0.619000 0.317000 + 19 ALA3 1 ALA DC DRUD -1.36530 0.400000 0 0.00000 0.00000 + 20 ALA3 1 ALA O OD2C1A 1.32050 15.5994 0 -0.579000 0.370000 + 21 ALA3 1 ALA DO DRUD -1.32050 0.400000 0 0.00000 0.00000 + 22 ALA3 1 ALA LPOA LPDO1 -0.311000 0.00000 -1 0.00000 0.00000 + 23 ALA3 1 ALA LPOB LPDO1 -0.227000 0.00000 -1 0.00000 0.00000 + 24 ALA3 1 ALA CB CD33A 2.23180 11.6110 0 -1.80400 2.08000 + 25 ALA3 1 ALA DCB DRUD -2.33080 0.400000 0 0.00000 0.00000 + 26 ALA3 1 ALA HB1 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 27 ALA3 1 ALA HB2 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 28 ALA3 1 ALA HB3 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 29 ALA3 2 ALA N ND2A2 1.95940 13.6070 0 -1.85800 0.126000 + 30 ALA3 2 ALA DN DRUD -2.36540 0.400000 0 0.00000 0.00000 + 31 ALA3 2 ALA HN HDP1A 0.296000 1.00800 0 0.00000 0.00000 + 32 ALA3 2 ALA CA CD31C 1.96170 11.6110 0 -1.06600 1.14100 + 33 ALA3 2 ALA DCA DRUD -1.79170 0.400000 0 0.00000 0.00000 + 34 ALA3 2 ALA HA HDA1A -0.170000E-01 1.00800 0 0.00000 0.00000 + 35 ALA3 2 ALA C CD2O1A 1.86030 11.6110 0 -0.619000 0.317000 + 36 ALA3 2 ALA DC DRUD -1.36530 0.400000 0 0.00000 0.00000 + 37 ALA3 2 ALA O OD2C1A 1.32050 15.5994 0 -0.579000 0.370000 + 38 ALA3 2 ALA DO DRUD -1.32050 0.400000 0 0.00000 0.00000 + 39 ALA3 2 ALA LPOA LPDO1 -0.311000 0.00000 -1 0.00000 0.00000 + 40 ALA3 2 ALA LPOB LPDO1 -0.227000 0.00000 -1 0.00000 0.00000 + 41 ALA3 2 ALA CB CD33A 2.23180 11.6110 0 -1.80400 2.08000 + 42 ALA3 2 ALA DCB DRUD -2.33080 0.400000 0 0.00000 0.00000 + 43 ALA3 2 ALA HB1 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 44 ALA3 2 ALA HB2 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 45 ALA3 2 ALA HB3 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 46 ALA3 3 ALA N ND2A2 1.95940 13.6070 0 -1.85800 0.126000 + 47 ALA3 3 ALA DN DRUD -2.36540 0.400000 0 0.00000 0.00000 + 48 ALA3 3 ALA HN HDP1A 0.296000 1.00800 0 0.00000 0.00000 + 49 ALA3 3 ALA CA CD31C 1.96170 11.6110 0 -1.06600 1.14100 + 50 ALA3 3 ALA DCA DRUD -1.79170 0.400000 0 0.00000 0.00000 + 51 ALA3 3 ALA HA HDA1A -0.170000E-01 1.00800 0 0.00000 0.00000 + 52 ALA3 3 ALA C CD2O1A 1.86030 11.6110 0 -0.619000 0.317000 + 53 ALA3 3 ALA DC DRUD -1.36530 0.400000 0 0.00000 0.00000 + 54 ALA3 3 ALA O OD2C1A 1.32050 15.5994 0 -0.579000 0.370000 + 55 ALA3 3 ALA DO DRUD -1.32050 0.400000 0 0.00000 0.00000 + 56 ALA3 3 ALA LPOA LPDO1 -0.311000 0.00000 -1 0.00000 0.00000 + 57 ALA3 3 ALA LPOB LPDO1 -0.227000 0.00000 -1 0.00000 0.00000 + 58 ALA3 3 ALA CB CD33A 2.23180 11.6110 0 -1.80400 2.08000 + 59 ALA3 3 ALA DCB DRUD -2.33080 0.400000 0 0.00000 0.00000 + 60 ALA3 3 ALA HB1 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 61 ALA3 3 ALA HB2 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 62 ALA3 3 ALA HB3 HDA3A 0.330000E-01 1.00800 0 0.00000 0.00000 + 63 ALA3 3 ALA NT ND2A2 2.03630 13.6070 0 -1.94200 0.250000 + 64 ALA3 3 ALA DNT DRUD -2.41830 0.400000 0 0.00000 0.00000 + 65 ALA3 3 ALA HNT HDP1A 0.272000 1.00800 0 0.00000 0.00000 + 66 ALA3 3 ALA CAT CD33C 2.16660 11.6110 0 -1.63900 2.12200 + 67 ALA3 3 ALA DCAT DRUD -2.22160 0.400000 0 0.00000 0.00000 + 68 ALA3 3 ALA HTC1 HDA3A 0.550000E-01 1.00800 0 0.00000 0.00000 + 69 ALA3 3 ALA HTC2 HDA3A 0.550000E-01 1.00800 0 0.00000 0.00000 + 70 ALA3 3 ALA HTC3 HDA3A 0.550000E-01 1.00800 0 0.00000 0.00000 + + 69 !NBOND: bonds + 1 2 6 7 8 9 12 6 + 1 6 6 8 1 3 1 4 + 1 5 8 10 8 11 12 13 + 15 16 18 19 20 21 24 25 + 12 15 15 18 18 29 15 17 + 12 14 18 20 15 24 24 26 + 24 27 24 28 20 22 20 23 + 29 30 32 33 35 36 37 38 + 41 42 29 32 32 35 35 46 + 32 34 29 31 35 37 32 41 + 41 43 41 44 41 45 37 39 + 37 40 46 47 49 50 52 53 + 54 55 58 59 46 49 49 52 + 49 51 46 48 52 54 49 58 + 58 60 58 61 58 62 54 56 + 54 57 63 64 66 67 52 63 + 63 65 63 66 66 68 66 69 + 66 70 + + 72 !NTHETA: angles + 3 1 4 3 1 5 3 1 6 + 4 1 5 4 1 6 5 1 6 + 1 6 8 1 6 12 8 6 12 + 6 12 14 6 12 15 14 12 15 + 12 15 17 12 15 18 12 15 24 + 17 15 18 17 15 24 18 15 24 + 15 18 20 15 18 29 20 18 29 + 15 24 26 15 24 27 15 24 28 + 26 24 27 26 24 28 27 24 28 + 18 29 31 18 29 32 31 29 32 + 29 32 34 29 32 35 29 32 41 + 34 32 35 34 32 41 35 32 41 + 32 35 37 32 35 46 37 35 46 + 32 41 43 32 41 44 32 41 45 + 43 41 44 43 41 45 44 41 45 + 35 46 48 35 46 49 48 46 49 + 46 49 51 46 49 52 46 49 58 + 51 49 52 51 49 58 52 49 58 + 49 52 54 49 52 63 54 52 63 + 49 58 60 49 58 61 49 58 62 + 60 58 61 60 58 62 61 58 62 + 52 63 65 52 63 66 65 63 66 + 63 66 68 63 66 69 63 66 70 + 68 66 69 68 66 70 69 66 70 + + 91 !NPHI: dihedrals + 1 6 12 14 1 6 12 15 + 3 1 6 8 3 1 6 12 + 4 1 6 8 4 1 6 12 + 5 1 6 8 5 1 6 12 + 6 12 15 17 6 12 15 18 + 6 12 15 24 8 6 12 14 + 8 6 12 15 12 15 18 20 + 12 15 18 29 12 15 24 26 + 12 15 24 27 12 15 24 28 + 14 12 15 17 14 12 15 18 + 14 12 15 24 15 18 29 31 + 15 18 29 32 17 15 18 20 + 17 15 18 29 17 15 24 26 + 17 15 24 27 17 15 24 28 + 18 15 24 26 18 15 24 27 + 18 15 24 28 18 29 32 34 + 18 29 32 35 18 29 32 41 + 20 18 15 24 20 18 29 31 + 20 18 29 32 24 15 18 29 + 29 32 35 37 29 32 35 46 + 29 32 41 43 29 32 41 44 + 29 32 41 45 31 29 32 34 + 31 29 32 35 31 29 32 41 + 32 35 46 48 32 35 46 49 + 34 32 35 37 34 32 35 46 + 34 32 41 43 34 32 41 44 + 34 32 41 45 35 32 41 43 + 35 32 41 44 35 32 41 45 + 35 46 49 51 35 46 49 52 + 35 46 49 58 37 35 32 41 + 37 35 46 48 37 35 46 49 + 41 32 35 46 46 49 52 54 + 46 49 52 63 46 49 58 60 + 46 49 58 61 46 49 58 62 + 48 46 49 51 48 46 49 52 + 48 46 49 58 49 52 63 65 + 49 52 63 66 51 49 52 54 + 51 49 52 63 51 49 58 60 + 51 49 58 61 51 49 58 62 + 52 49 58 60 52 49 58 61 + 52 49 58 62 52 63 66 68 + 52 63 66 69 52 63 66 70 + 54 52 49 58 54 52 63 65 + 54 52 63 66 58 49 52 63 + 65 63 66 68 65 63 66 69 + 65 63 66 70 + + 8 !NIMPHI: impropers + 6 1 12 8 12 6 15 14 + 18 15 29 20 35 32 46 37 + 29 18 32 31 46 35 49 48 + 52 49 63 54 63 52 66 65 + + 4 !NDON: donors + 12 14 29 31 46 48 63 65 + + 4 !NACC: acceptors + 8 6 20 18 37 35 54 52 + + 0 !NNB + + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 + + 7 0 !NGRP NST2 + 0 1 0 11 1 0 23 1 0 + 28 1 0 40 1 0 45 1 0 + 57 1 0 + + 8 32 !NUMLP NUMLPH + 3 1 F 0.300000 91.0000 0.00000 + 3 5 F 0.300000 91.0000 180.000 + 3 9 F 0.300000 91.0000 0.00000 + 3 13 F 0.300000 91.0000 180.000 + 3 17 F 0.300000 91.0000 0.00000 + 3 21 F 0.300000 91.0000 180.000 + 3 25 F 0.300000 91.0000 0.00000 + 3 29 F 0.300000 91.0000 180.000 + 10 8 6 1 11 8 6 1 + 22 20 18 15 23 20 18 15 + 39 37 35 32 40 37 35 32 + 56 54 52 49 57 54 52 49 + + 4 !NUMANISO + 123.559 -46.4888 -16.1883 + 123.559 -46.4888 -16.1883 + 123.559 -46.4888 -16.1883 + 123.559 -46.4888 -16.1883 + 8 6 10 11 20 18 22 23 + 37 35 39 40 54 52 56 57 + + 3 !NCRTERM: cross-terms + 6 12 15 18 12 15 18 29 + 18 29 32 35 29 32 35 46 + 35 46 49 52 46 49 52 63 + + 0 !AUTOGEN + 0 15 0 0 0 0 15 0 + 15 15 15 0 15 0 0 15 + 0 0 15 0 15 15 15 0 + 15 0 0 0 0 15 0 0 + 15 0 0 15 0 15 15 15 + 0 15 0 0 0 0 15 0 + 0 15 0 0 15 0 15 15 + 15 0 15 0 0 0 0 15 + 0 0 15 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0