diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 8c073b0ad..7651113f4 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -144,19 +144,18 @@ set(_FILES orbital orbits orient orthog output overlap params paths pbstuf pdb phipsi picalc piorbs pistuf pitors pme pmestuf pmpb polar polgrp polopt polpcg polpot poltcg polymer potent potfit predict pressure - prmkey promo prtarc prtcif prtdyn prterr prtfrc prtint prtmol2 prtpdb - prtprm prtseq prtuind prtvel prtxyz ptable qmstuf qrsolve quatfit - random rattle readcart readcif readdcd readdyn readgau readgdma - readint readmbis readmol readmol2 readpdb readprm readseq readxyz - refer repel replica reppot resdue respa restrn rgddyn rgdstep richmond - rigid ring rings rmsfit rotbnd rotlist rotpole rxnfld rxnpot scales - sdstep search sequen server setprm settle shake shapes shunt sigmoid - simplex sizes sktstuf socket solpot solute sort square stodyn strbnd - strtor suffix surface surfatom switch syntrn tarray tcgstuf temper - tettor tettors titles tncg torphase torpot torque tors torsions tortor - tree trimtext tritor tritors unionball unitcell units uprior urey - urypot usage valfit vdw vdwpot verlet version vibs virial volume warp - xtals xyzatm zatom zclose zcoord + prmkey promo prtarc prtcif prtdyn prterr prtint prtmol2 prtpdb prtprm + prtseq prtxyz ptable qmstuf qrsolve quatfit random rattle readcart + readcif readdcd readdyn readgau readgdma readint readmbis readmol + readmol2 readpdb readprm readseq readxyz refer repel replica reppot + resdue respa restrn rgddyn rgdstep richmond rigid ring rings rmsfit + rotbnd rotlist rotpole rxnfld rxnpot scales sdstep search sequen server + setprm settle shake shapes shunt sigmoid simplex sizes sktstuf socket + solpot solute sort square stodyn strbnd strtor suffix surface surfatom + switch syntrn tarray tcgstuf temper tettor tettors titles tncg + torphase torpot torque tors torsions tortor tree trimtext tritor tritors + unionball unitcell units uprior urey urypot usage valfit vdw vdwpot + verlet version vibs virial volume warp xtals xyzatm zatom zclose zcoord ) set(_LIB "") diff --git a/interface/CMakeLists.txt b/interface/CMakeLists.txt index c72133e3d..edb55fc97 100644 --- a/interface/CMakeLists.txt +++ b/interface/CMakeLists.txt @@ -166,6 +166,7 @@ add_library (tinkerObjF OBJECT ../source/tortor.f ../source/tree.f ../source/tritor.f +../source/uatom.f ../source/units.f ../source/uprior.f ../source/urey.f @@ -482,14 +483,11 @@ add_library (tinkerObjF OBJECT ../source/prtcif.f ../source/prtdyn.f ../source/prterr.f -../source/prtfrc.f ../source/prtint.f ../source/prtmol2.f ../source/prtpdb.f ../source/prtprm.f ../source/prtseq.f -../source/prtuind.f -../source/prtvel.f ../source/prtxyz.f ../source/qrsolve.f ../source/quatfit.f @@ -542,6 +540,7 @@ add_library (tinkerObjF OBJECT ../source/trimtext.f ../source/tritors.f ../source/unionball.f +../source/uniquetyp.f ../source/unitcell.f ../source/verlet.f ../source/version.f diff --git a/interface/c/tinker/detail/boxes.hh b/interface/c/tinker/detail/boxes.hh index b5de41ea5..b59e1acde 100644 --- a/interface/c/tinker/detail/boxes.hh +++ b/interface/c/tinker/detail/boxes.hh @@ -24,6 +24,9 @@ extern double TINKER_MOD(boxes, gamma_sin); extern double TINKER_MOD(boxes, gamma_cos); extern double TINKER_MOD(boxes, beta_term); extern double TINKER_MOD(boxes, gamma_term); +extern double TINKER_MOD(boxes, xcenter); +extern double TINKER_MOD(boxes, ycenter); +extern double TINKER_MOD(boxes, zcenter); extern double TINKER_MOD(boxes, lvec)[3][3]; extern double TINKER_MOD(boxes, recip)[3][3]; extern int TINKER_MOD(boxes, orthogonal); diff --git a/interface/c/tinker/detail/extfld.hh b/interface/c/tinker/detail/extfld.hh index 6e431d2d0..b49d6e581 100644 --- a/interface/c/tinker/detail/extfld.hh +++ b/interface/c/tinker/detail/extfld.hh @@ -5,8 +5,11 @@ #ifdef __cplusplus extern "C" { #endif +extern double TINKER_MOD(extfld, exfreq); extern double TINKER_MOD(extfld, exfld)[3]; +extern double TINKER_MOD(extfld, texfld)[3]; extern int TINKER_MOD(extfld, use_exfld); +extern int TINKER_MOD(extfld, use_exfreq); #ifdef __cplusplus } #endif diff --git a/interface/c/tinker/detail/moment.hh b/interface/c/tinker/detail/moment.hh index 993c087de..a3bbcbfe0 100644 --- a/interface/c/tinker/detail/moment.hh +++ b/interface/c/tinker/detail/moment.hh @@ -20,6 +20,7 @@ extern double TINKER_MOD(moment, yzqpl); extern double TINKER_MOD(moment, zxqpl); extern double TINKER_MOD(moment, zyqpl); extern double TINKER_MOD(moment, zzqpl); +extern int* TINKER_MOD(moment, momuse); #ifdef __cplusplus } #endif diff --git a/interface/c/tinker/detail/output.hh b/interface/c/tinker/detail/output.hh index fdf26b2f0..f526c69bc 100644 --- a/interface/c/tinker/detail/output.hh +++ b/interface/c/tinker/detail/output.hh @@ -5,16 +5,30 @@ #ifdef __cplusplus extern "C" { #endif +extern int TINKER_MOD(output, nonly); +extern int* TINKER_MOD(output, ionly); +extern int* TINKER_MOD(output, ionlyinv); +extern double* TINKER_MOD(output, print3n); extern int TINKER_MOD(output, archive); extern int TINKER_MOD(output, binary); extern int TINKER_MOD(output, noversion); extern int TINKER_MOD(output, overwrite); +extern int TINKER_MOD(output, coordsave); +extern int TINKER_MOD(output, dynsave); extern int TINKER_MOD(output, cyclesave); +extern int TINKER_MOD(output, onlysave); extern int TINKER_MOD(output, arcsave); extern int TINKER_MOD(output, dcdsave); extern int TINKER_MOD(output, velsave); extern int TINKER_MOD(output, frcsave); extern int TINKER_MOD(output, uindsave); +extern int TINKER_MOD(output, ustcsave); +extern int TINKER_MOD(output, uchgsave); +extern int TINKER_MOD(output, usyssave); +extern int TINKER_MOD(output, vsyssave); +extern int TINKER_MOD(output, udirsave); +extern int TINKER_MOD(output, defsave); +extern int TINKER_MOD(output, tefsave); extern char TINKER_MOD(output, coordtype)[9]; #ifdef __cplusplus } diff --git a/interface/c/tinker/detail/uatom.hh b/interface/c/tinker/detail/uatom.hh new file mode 100644 index 000000000..40d76e44e --- /dev/null +++ b/interface/c/tinker/detail/uatom.hh @@ -0,0 +1,17 @@ +#pragma once + +#include "macro.hh" +#include "sizes.hh" + +#ifdef __cplusplus +extern "C" { +#endif +extern int TINKER_MOD(uatom, nunique); +extern int TINKER_MOD(uatom, utype)[TINKER_MOD__maxtyp]; +extern int TINKER_MOD(uatom, utypeinv)[TINKER_MOD__maxtyp]; +extern double TINKER_MOD(uatom, utv1)[TINKER_MOD__maxtyp][3]; +extern double TINKER_MOD(uatom, utv2)[TINKER_MOD__maxtyp][3]; +extern double TINKER_MOD(uatom, utv3)[TINKER_MOD__maxtyp][3]; +#ifdef __cplusplus +} +#endif diff --git a/interface/c/tinker/modc.h b/interface/c/tinker/modc.h index b924d39b3..4bb8dea0b 100644 --- a/interface/c/tinker/modc.h +++ b/interface/c/tinker/modc.h @@ -162,6 +162,7 @@ #include "detail/tortor.hh" #include "detail/tree.hh" #include "detail/tritor.hh" +#include "detail/uatom.hh" #include "detail/units.hh" #include "detail/uprior.hh" #include "detail/urey.hh" diff --git a/interface/c/tinker/routines.h b/interface/c/tinker/routines.h index 4d57f4926..c96715cdc 100644 --- a/interface/c/tinker/routines.h +++ b/interface/c/tinker/routines.h @@ -20,6 +20,10 @@ void tinkerFortranRuntimeEnd(); // active.f void active_(); #define tinker_f_active active_ +void saveonly_(); +#define tinker_f_saveonly saveonly_ +void msystem_(); +#define tinker_f_msystem msystem_ // alterchg.f void alterchg_(); @@ -108,6 +112,8 @@ void calendar_(int* year, int* month, int* day, int* hour, int* minute, int* sec // center.f void center_(int* n1, double* x1, double* y1, double* z1, int* n2, double* x2, double* y2, double* z2, double* xmid, double* ymid, double* zmid); #define tinker_f_center center_ +void compcent_(double* xcm, double* ycm, double* zcm); +#define tinker_f_compcent compcent_ // chkpole.f void chkpole_(); @@ -1970,6 +1976,8 @@ void moments_(char* mode, tinker_fchar_len_t mode_cap); inline void tinker_f_moments(tinker_fchars mode) { return moments_(mode.string, mode.capacity); } +void dmoments_(double* xm, double* ym, double* zm, double* xustc, double* yustc, double* zustc, double* xuind, double* yuind, double* zuind, double* xuchg, double* yuchg, double* zuchg); +#define tinker_f_dmoments dmoments_ // mutate.f void mutate_(); @@ -2247,12 +2255,6 @@ void prtdyn_(); void prterr_(); #define tinker_f_prterr prterr_ -// prtfrc.f -void prtfrc_(int* ifrc); -#define tinker_f_prtfrc prtfrc_ -void prtdcdf_(int* idcd, int* first); -#define tinker_f_prtdcdf prtdcdf_ - // prtint.f void prtint_(int* izmt); #define tinker_f_prtint prtint_ @@ -2277,23 +2279,23 @@ void prtprm_(int* itxt); void prtseq_(int* iseq); #define tinker_f_prtseq prtseq_ -// prtuind.f -void prtuind_(int* iind); -#define tinker_f_prtuind prtuind_ -void prtdcdu_(int* idcd, int* first); -#define tinker_f_prtdcdu prtdcdu_ - -// prtvel.f -void prtvel_(int* ivel); -#define tinker_f_prtvel prtvel_ -void prtdcdv_(int* idcd, int* first); -#define tinker_f_prtdcdv prtdcdv_ - // prtxyz.f void prtxyz_(int* ixyz); #define tinker_f_prtxyz prtxyz_ +void prtvec3_(int* iunit, char* mode, tinker_fchar_len_t mode_cap); +inline void tinker_f_prtvec3(int* iunit, tinker_fchars mode) { + return prtvec3_(iunit, mode.string, mode.capacity); +} +void copyvec3_(double* print3n, char* mode, tinker_fchar_len_t mode_cap); +inline void tinker_f_copyvec3(double* print3n, tinker_fchars mode) { + return copyvec3_(print3n, mode.string, mode.capacity); +} void prtdcd_(int* idcd, int* first); #define tinker_f_prtdcd prtdcd_ +void prtdcdv3_(int* idcd, int* first, char* mode, tinker_fchar_len_t mode_cap); +inline void tinker_f_prtdcdv3(int* idcd, int* first, tinker_fchars mode) { + return prtdcdv3_(idcd, first, mode.string, mode.capacity); +} // qrsolve.f void qrfact_(int* n, int* m, double* a, int* pivot, int* ipvt, double* rdiag); @@ -2913,6 +2915,12 @@ void addbogus_(double* bcoord, double* brad); void tetra_volume_(double* r12sq, double* r13sq, double* r14sq, double* r23sq, double* r24sq, double* r34sq, double* vol); #define tinker_f_tetra_volume tetra_volume_ +// uniquetyp.f +void uniquetyp_(); +#define tinker_f_uniquetyp uniquetyp_ +void velunique_(); +#define tinker_f_velunique velunique_ + // unitcell.f void unitcell_(); #define tinker_f_unitcell unitcell_ diff --git a/interface/cpp/tinker/detail/boxes.hh b/interface/cpp/tinker/detail/boxes.hh index 23b5283fc..ad013fbae 100644 --- a/interface/cpp/tinker/detail/boxes.hh +++ b/interface/cpp/tinker/detail/boxes.hh @@ -22,6 +22,9 @@ extern double& gamma_sin; extern double& gamma_cos; extern double& beta_term; extern double& gamma_term; +extern double& xcenter; +extern double& ycenter; +extern double& zcenter; extern double (&lvec)[3][3]; extern double (&recip)[3][3]; extern int& orthogonal; @@ -53,6 +56,9 @@ extern "C" double TINKER_MOD(boxes, gamma_sin); extern "C" double TINKER_MOD(boxes, gamma_cos); extern "C" double TINKER_MOD(boxes, beta_term); extern "C" double TINKER_MOD(boxes, gamma_term); +extern "C" double TINKER_MOD(boxes, xcenter); +extern "C" double TINKER_MOD(boxes, ycenter); +extern "C" double TINKER_MOD(boxes, zcenter); extern "C" double TINKER_MOD(boxes, lvec)[3][3]; extern "C" double TINKER_MOD(boxes, recip)[3][3]; extern "C" int TINKER_MOD(boxes, orthogonal); @@ -83,6 +89,9 @@ double& gamma_sin = TINKER_MOD(boxes, gamma_sin); double& gamma_cos = TINKER_MOD(boxes, gamma_cos); double& beta_term = TINKER_MOD(boxes, beta_term); double& gamma_term = TINKER_MOD(boxes, gamma_term); +double& xcenter = TINKER_MOD(boxes, xcenter); +double& ycenter = TINKER_MOD(boxes, ycenter); +double& zcenter = TINKER_MOD(boxes, zcenter); double (&lvec)[3][3] = TINKER_MOD(boxes, lvec); double (&recip)[3][3] = TINKER_MOD(boxes, recip); int& orthogonal = TINKER_MOD(boxes, orthogonal); diff --git a/interface/cpp/tinker/detail/extfld.hh b/interface/cpp/tinker/detail/extfld.hh index fe6555c9f..dc2e8e47b 100644 --- a/interface/cpp/tinker/detail/extfld.hh +++ b/interface/cpp/tinker/detail/extfld.hh @@ -3,14 +3,23 @@ #include "macro.hh" namespace tinker { namespace extfld { +extern double& exfreq; extern double (&exfld)[3]; +extern double (&texfld)[3]; extern int& use_exfld; +extern int& use_exfreq; #ifdef TINKER_FORTRAN_MODULE_CPP +extern "C" double TINKER_MOD(extfld, exfreq); extern "C" double TINKER_MOD(extfld, exfld)[3]; +extern "C" double TINKER_MOD(extfld, texfld)[3]; extern "C" int TINKER_MOD(extfld, use_exfld); +extern "C" int TINKER_MOD(extfld, use_exfreq); +double& exfreq = TINKER_MOD(extfld, exfreq); double (&exfld)[3] = TINKER_MOD(extfld, exfld); +double (&texfld)[3] = TINKER_MOD(extfld, texfld); int& use_exfld = TINKER_MOD(extfld, use_exfld); +int& use_exfreq = TINKER_MOD(extfld, use_exfreq); #endif } } diff --git a/interface/cpp/tinker/detail/moment.hh b/interface/cpp/tinker/detail/moment.hh index 8c4d3bff6..a7483cab0 100644 --- a/interface/cpp/tinker/detail/moment.hh +++ b/interface/cpp/tinker/detail/moment.hh @@ -18,6 +18,7 @@ extern double& yzqpl; extern double& zxqpl; extern double& zyqpl; extern double& zzqpl; +extern int*& momuse; #ifdef TINKER_FORTRAN_MODULE_CPP extern "C" double TINKER_MOD(moment, netchg); @@ -35,6 +36,7 @@ extern "C" double TINKER_MOD(moment, yzqpl); extern "C" double TINKER_MOD(moment, zxqpl); extern "C" double TINKER_MOD(moment, zyqpl); extern "C" double TINKER_MOD(moment, zzqpl); +extern "C" int* TINKER_MOD(moment, momuse); double& netchg = TINKER_MOD(moment, netchg); double& netdpl = TINKER_MOD(moment, netdpl); @@ -51,5 +53,6 @@ double& yzqpl = TINKER_MOD(moment, yzqpl); double& zxqpl = TINKER_MOD(moment, zxqpl); double& zyqpl = TINKER_MOD(moment, zyqpl); double& zzqpl = TINKER_MOD(moment, zzqpl); +int*& momuse = TINKER_MOD(moment, momuse); #endif } } diff --git a/interface/cpp/tinker/detail/output.hh b/interface/cpp/tinker/detail/output.hh index a71ad373c..4b962e389 100644 --- a/interface/cpp/tinker/detail/output.hh +++ b/interface/cpp/tinker/detail/output.hh @@ -3,41 +3,83 @@ #include "macro.hh" namespace tinker { namespace output { +extern int& nonly; +extern int*& ionly; +extern int*& ionlyinv; +extern double*& print3n; extern int& archive; extern int& binary; extern int& noversion; extern int& overwrite; +extern int& coordsave; +extern int& dynsave; extern int& cyclesave; +extern int& onlysave; extern int& arcsave; extern int& dcdsave; extern int& velsave; extern int& frcsave; extern int& uindsave; +extern int& ustcsave; +extern int& uchgsave; +extern int& usyssave; +extern int& vsyssave; +extern int& udirsave; +extern int& defsave; +extern int& tefsave; extern char (&coordtype)[9]; #ifdef TINKER_FORTRAN_MODULE_CPP +extern "C" int TINKER_MOD(output, nonly); +extern "C" int* TINKER_MOD(output, ionly); +extern "C" int* TINKER_MOD(output, ionlyinv); +extern "C" double* TINKER_MOD(output, print3n); extern "C" int TINKER_MOD(output, archive); extern "C" int TINKER_MOD(output, binary); extern "C" int TINKER_MOD(output, noversion); extern "C" int TINKER_MOD(output, overwrite); +extern "C" int TINKER_MOD(output, coordsave); +extern "C" int TINKER_MOD(output, dynsave); extern "C" int TINKER_MOD(output, cyclesave); +extern "C" int TINKER_MOD(output, onlysave); extern "C" int TINKER_MOD(output, arcsave); extern "C" int TINKER_MOD(output, dcdsave); extern "C" int TINKER_MOD(output, velsave); extern "C" int TINKER_MOD(output, frcsave); extern "C" int TINKER_MOD(output, uindsave); +extern "C" int TINKER_MOD(output, ustcsave); +extern "C" int TINKER_MOD(output, uchgsave); +extern "C" int TINKER_MOD(output, usyssave); +extern "C" int TINKER_MOD(output, vsyssave); +extern "C" int TINKER_MOD(output, udirsave); +extern "C" int TINKER_MOD(output, defsave); +extern "C" int TINKER_MOD(output, tefsave); extern "C" char TINKER_MOD(output, coordtype)[9]; +int& nonly = TINKER_MOD(output, nonly); +int*& ionly = TINKER_MOD(output, ionly); +int*& ionlyinv = TINKER_MOD(output, ionlyinv); +double*& print3n = TINKER_MOD(output, print3n); int& archive = TINKER_MOD(output, archive); int& binary = TINKER_MOD(output, binary); int& noversion = TINKER_MOD(output, noversion); int& overwrite = TINKER_MOD(output, overwrite); +int& coordsave = TINKER_MOD(output, coordsave); +int& dynsave = TINKER_MOD(output, dynsave); int& cyclesave = TINKER_MOD(output, cyclesave); +int& onlysave = TINKER_MOD(output, onlysave); int& arcsave = TINKER_MOD(output, arcsave); int& dcdsave = TINKER_MOD(output, dcdsave); int& velsave = TINKER_MOD(output, velsave); int& frcsave = TINKER_MOD(output, frcsave); int& uindsave = TINKER_MOD(output, uindsave); +int& ustcsave = TINKER_MOD(output, ustcsave); +int& uchgsave = TINKER_MOD(output, uchgsave); +int& usyssave = TINKER_MOD(output, usyssave); +int& vsyssave = TINKER_MOD(output, vsyssave); +int& udirsave = TINKER_MOD(output, udirsave); +int& defsave = TINKER_MOD(output, defsave); +int& tefsave = TINKER_MOD(output, tefsave); char (&coordtype)[9] = TINKER_MOD(output, coordtype); #endif } } diff --git a/interface/cpp/tinker/detail/uatom.hh b/interface/cpp/tinker/detail/uatom.hh new file mode 100644 index 000000000..f6d853c43 --- /dev/null +++ b/interface/cpp/tinker/detail/uatom.hh @@ -0,0 +1,31 @@ +#pragma once + +#include "macro.hh" +#include "sizes.hh" + +namespace tinker { namespace uatom { +using namespace sizes; + +extern int& nunique; +extern int (&utype)[maxtyp]; +extern int (&utypeinv)[maxtyp]; +extern double (&utv1)[maxtyp][3]; +extern double (&utv2)[maxtyp][3]; +extern double (&utv3)[maxtyp][3]; + +#ifdef TINKER_FORTRAN_MODULE_CPP +extern "C" int TINKER_MOD(uatom, nunique); +extern "C" int TINKER_MOD(uatom, utype)[maxtyp]; +extern "C" int TINKER_MOD(uatom, utypeinv)[maxtyp]; +extern "C" double TINKER_MOD(uatom, utv1)[maxtyp][3]; +extern "C" double TINKER_MOD(uatom, utv2)[maxtyp][3]; +extern "C" double TINKER_MOD(uatom, utv3)[maxtyp][3]; + +int& nunique = TINKER_MOD(uatom, nunique); +int (&utype)[maxtyp] = TINKER_MOD(uatom, utype); +int (&utypeinv)[maxtyp] = TINKER_MOD(uatom, utypeinv); +double (&utv1)[maxtyp][3] = TINKER_MOD(uatom, utv1); +double (&utv2)[maxtyp][3] = TINKER_MOD(uatom, utv2); +double (&utv3)[maxtyp][3] = TINKER_MOD(uatom, utv3); +#endif +} } diff --git a/interface/cpp/tinker/modcpp.h b/interface/cpp/tinker/modcpp.h index b924d39b3..4bb8dea0b 100644 --- a/interface/cpp/tinker/modcpp.h +++ b/interface/cpp/tinker/modcpp.h @@ -162,6 +162,7 @@ #include "detail/tortor.hh" #include "detail/tree.hh" #include "detail/tritor.hh" +#include "detail/uatom.hh" #include "detail/units.hh" #include "detail/uprior.hh" #include "detail/urey.hh" diff --git a/interface/cpp/tinker/routines.h b/interface/cpp/tinker/routines.h index 4d57f4926..c96715cdc 100644 --- a/interface/cpp/tinker/routines.h +++ b/interface/cpp/tinker/routines.h @@ -20,6 +20,10 @@ void tinkerFortranRuntimeEnd(); // active.f void active_(); #define tinker_f_active active_ +void saveonly_(); +#define tinker_f_saveonly saveonly_ +void msystem_(); +#define tinker_f_msystem msystem_ // alterchg.f void alterchg_(); @@ -108,6 +112,8 @@ void calendar_(int* year, int* month, int* day, int* hour, int* minute, int* sec // center.f void center_(int* n1, double* x1, double* y1, double* z1, int* n2, double* x2, double* y2, double* z2, double* xmid, double* ymid, double* zmid); #define tinker_f_center center_ +void compcent_(double* xcm, double* ycm, double* zcm); +#define tinker_f_compcent compcent_ // chkpole.f void chkpole_(); @@ -1970,6 +1976,8 @@ void moments_(char* mode, tinker_fchar_len_t mode_cap); inline void tinker_f_moments(tinker_fchars mode) { return moments_(mode.string, mode.capacity); } +void dmoments_(double* xm, double* ym, double* zm, double* xustc, double* yustc, double* zustc, double* xuind, double* yuind, double* zuind, double* xuchg, double* yuchg, double* zuchg); +#define tinker_f_dmoments dmoments_ // mutate.f void mutate_(); @@ -2247,12 +2255,6 @@ void prtdyn_(); void prterr_(); #define tinker_f_prterr prterr_ -// prtfrc.f -void prtfrc_(int* ifrc); -#define tinker_f_prtfrc prtfrc_ -void prtdcdf_(int* idcd, int* first); -#define tinker_f_prtdcdf prtdcdf_ - // prtint.f void prtint_(int* izmt); #define tinker_f_prtint prtint_ @@ -2277,23 +2279,23 @@ void prtprm_(int* itxt); void prtseq_(int* iseq); #define tinker_f_prtseq prtseq_ -// prtuind.f -void prtuind_(int* iind); -#define tinker_f_prtuind prtuind_ -void prtdcdu_(int* idcd, int* first); -#define tinker_f_prtdcdu prtdcdu_ - -// prtvel.f -void prtvel_(int* ivel); -#define tinker_f_prtvel prtvel_ -void prtdcdv_(int* idcd, int* first); -#define tinker_f_prtdcdv prtdcdv_ - // prtxyz.f void prtxyz_(int* ixyz); #define tinker_f_prtxyz prtxyz_ +void prtvec3_(int* iunit, char* mode, tinker_fchar_len_t mode_cap); +inline void tinker_f_prtvec3(int* iunit, tinker_fchars mode) { + return prtvec3_(iunit, mode.string, mode.capacity); +} +void copyvec3_(double* print3n, char* mode, tinker_fchar_len_t mode_cap); +inline void tinker_f_copyvec3(double* print3n, tinker_fchars mode) { + return copyvec3_(print3n, mode.string, mode.capacity); +} void prtdcd_(int* idcd, int* first); #define tinker_f_prtdcd prtdcd_ +void prtdcdv3_(int* idcd, int* first, char* mode, tinker_fchar_len_t mode_cap); +inline void tinker_f_prtdcdv3(int* idcd, int* first, tinker_fchars mode) { + return prtdcdv3_(idcd, first, mode.string, mode.capacity); +} // qrsolve.f void qrfact_(int* n, int* m, double* a, int* pivot, int* ipvt, double* rdiag); @@ -2913,6 +2915,12 @@ void addbogus_(double* bcoord, double* brad); void tetra_volume_(double* r12sq, double* r13sq, double* r14sq, double* r23sq, double* r24sq, double* r34sq, double* vol); #define tinker_f_tetra_volume tetra_volume_ +// uniquetyp.f +void uniquetyp_(); +#define tinker_f_uniquetyp uniquetyp_ +void velunique_(); +#define tinker_f_velunique velunique_ + // unitcell.f void unitcell_(); #define tinker_f_unitcell unitcell_ diff --git a/linux/gfortran/compile.make b/linux/gfortran/compile.make index 43b7be2d2..2f25a8602 100755 --- a/linux/gfortran/compile.make +++ b/linux/gfortran/compile.make @@ -528,14 +528,11 @@ gfortran -c -O3 -fopenmp prtarc.f gfortran -c -O3 -fopenmp prtcif.f gfortran -c -O3 -fopenmp prtdyn.f gfortran -c -O3 -fopenmp prterr.f -gfortran -c -O3 -fopenmp prtfrc.f gfortran -c -O3 -fopenmp prtint.f gfortran -c -O3 -fopenmp prtmol2.f gfortran -c -O3 -fopenmp prtpdb.f gfortran -c -O3 -fopenmp prtprm.f gfortran -c -O3 -fopenmp prtseq.f -gfortran -c -O3 -fopenmp prtuind.f -gfortran -c -O3 -fopenmp prtvel.f gfortran -c -O3 -fopenmp prtxyz.f gfortran -c -O3 -fopenmp pss.f gfortran -c -O3 -fopenmp pssrigid.f diff --git a/linux/gfortran/debug.make b/linux/gfortran/debug.make index 229c25dc2..5687098ba 100755 --- a/linux/gfortran/debug.make +++ b/linux/gfortran/debug.make @@ -528,14 +528,11 @@ gfortran -c -Og -g -Wall prtarc.f gfortran -c -Og -g -Wall prtcif.f gfortran -c -Og -g -Wall prtdyn.f gfortran -c -Og -g -Wall prterr.f -gfortran -c -Og -g -Wall prtfrc.f gfortran -c -Og -g -Wall prtint.f gfortran -c -Og -g -Wall prtmol2.f gfortran -c -Og -g -Wall prtpdb.f gfortran -c -Og -g -Wall prtprm.f gfortran -c -Og -g -Wall prtseq.f -gfortran -c -Og -g -Wall prtuind.f -gfortran -c -Og -g -Wall prtvel.f gfortran -c -Og -g -Wall prtxyz.f gfortran -c -Og -g -Wall pss.f gfortran -c -Og -g -Wall pssrigid.f diff --git a/linux/gfortran/generic.make b/linux/gfortran/generic.make index 73e0a4e83..2334d24bc 100755 --- a/linux/gfortran/generic.make +++ b/linux/gfortran/generic.make @@ -528,14 +528,11 @@ gfortran -c -O3 -mtune=generic -fopenmp prtarc.f gfortran -c -O3 -mtune=generic -fopenmp prtcif.f gfortran -c -O3 -mtune=generic -fopenmp prtdyn.f gfortran -c -O3 -mtune=generic -fopenmp prterr.f -gfortran -c -O3 -mtune=generic -fopenmp prtfrc.f gfortran -c -O3 -mtune=generic -fopenmp prtint.f gfortran -c -O3 -mtune=generic -fopenmp prtmol2.f gfortran -c -O3 -mtune=generic -fopenmp prtpdb.f gfortran -c -O3 -mtune=generic -fopenmp prtprm.f gfortran -c -O3 -mtune=generic -fopenmp prtseq.f -gfortran -c -O3 -mtune=generic -fopenmp prtuind.f -gfortran -c -O3 -mtune=generic -fopenmp prtvel.f gfortran -c -O3 -mtune=generic -fopenmp prtxyz.f gfortran -c -O3 -mtune=generic -fopenmp pss.f gfortran -c -O3 -mtune=generic -fopenmp pssrigid.f diff --git a/linux/gfortran/library.make b/linux/gfortran/library.make index 73f4ba780..1675a4ad0 100755 --- a/linux/gfortran/library.make +++ b/linux/gfortran/library.make @@ -439,14 +439,11 @@ prtarc.o \ prtcif.o \ prtdyn.o \ prterr.o \ -prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ -prtuind.o \ -prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ diff --git a/linux/intel/compile.make b/linux/intel/compile.make index a09d58720..0bdf3eacf 100755 --- a/linux/intel/compile.make +++ b/linux/intel/compile.make @@ -528,14 +528,11 @@ ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtarc.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtcif.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtdyn.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prterr.f -ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtfrc.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtint.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtmol2.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtpdb.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtprm.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtseq.f -ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtuind.f -ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtvel.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp prtxyz.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp pss.f ifort -c -O3 -no-ipo -no-prec-div -qopenmp pssrigid.f diff --git a/linux/intel/compserial.make b/linux/intel/compserial.make index 65e668ffc..a246f5063 100755 --- a/linux/intel/compserial.make +++ b/linux/intel/compserial.make @@ -528,14 +528,11 @@ ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtarc.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtcif.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtdyn.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prterr.f -ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtfrc.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtint.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtmol2.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtpdb.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtprm.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtseq.f -ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtuind.f -ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtvel.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 prtxyz.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 pss.f ifort -c -O3 -no-ipo -no-prec-div -vec-report0 pssrigid.f diff --git a/linux/intel/debug.make b/linux/intel/debug.make index 06cae9115..986e44199 100755 --- a/linux/intel/debug.make +++ b/linux/intel/debug.make @@ -528,14 +528,11 @@ ifort -c -g -warn all -check all prtarc.f ifort -c -g -warn all -check all prtcif.f ifort -c -g -warn all -check all prtdyn.f ifort -c -g -warn all -check all prterr.f -ifort -c -g -warn all -check all prtfrc.f ifort -c -g -warn all -check all prtint.f ifort -c -g -warn all -check all prtmol2.f ifort -c -g -warn all -check all prtpdb.f ifort -c -g -warn all -check all prtprm.f ifort -c -g -warn all -check all prtseq.f -ifort -c -g -warn all -check all prtuind.f -ifort -c -g -warn all -check all prtvel.f ifort -c -g -warn all -check all prtxyz.f ifort -c -g -warn all -check all pss.f ifort -c -g -warn all -check all pssrigid.f diff --git a/linux/intel/generic.make b/linux/intel/generic.make index 70a842696..374d81da5 100755 --- a/linux/intel/generic.make +++ b/linux/intel/generic.make @@ -528,14 +528,11 @@ ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtarc.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtcif.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtdyn.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prterr.f -ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtfrc.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtint.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtmol2.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtpdb.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtprm.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtseq.f -ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtuind.f -ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtvel.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp prtxyz.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp pss.f ifort -c -O3 -msse3 -no-ipo -no-prec-div -qopenmp pssrigid.f diff --git a/linux/intel/library.make b/linux/intel/library.make index 47c0f8549..1f3e89849 100755 --- a/linux/intel/library.make +++ b/linux/intel/library.make @@ -439,14 +439,11 @@ prtarc.o \ prtcif.o \ prtdyn.o \ prterr.o \ -prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ -prtuind.o \ -prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ diff --git a/macos/flang/compile.make b/macos/flang/compile.make index b75a96c4b..ee7fa5bdc 100755 --- a/macos/flang/compile.make +++ b/macos/flang/compile.make @@ -531,14 +531,11 @@ flang -c -O3 -fopenmp prtarc.f flang -c -O3 -fopenmp prtcif.f flang -c -O3 -fopenmp prtdyn.f flang -c -O3 -fopenmp prterr.f -flang -c -O3 -fopenmp prtfrc.f flang -c -O3 -fopenmp prtint.f flang -c -O3 -fopenmp prtmol2.f flang -c -O3 -fopenmp prtpdb.f flang -c -O3 -fopenmp prtprm.f flang -c -O3 -fopenmp prtseq.f -flang -c -O3 -fopenmp prtuind.f -flang -c -O3 -fopenmp prtvel.f flang -c -O3 -fopenmp prtxyz.f flang -c -O3 -fopenmp pss.f flang -c -O3 -fopenmp pssrigid.f diff --git a/macos/flang/library.make b/macos/flang/library.make index 0d75d2813..6891a8eec 100755 --- a/macos/flang/library.make +++ b/macos/flang/library.make @@ -439,14 +439,11 @@ prtarc.o \ prtcif.o \ prtdyn.o \ prterr.o \ -prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ -prtuind.o \ -prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ diff --git a/macos/gfortran/compile.make b/macos/gfortran/compile.make index ced111e26..36f0a2da7 100755 --- a/macos/gfortran/compile.make +++ b/macos/gfortran/compile.make @@ -528,14 +528,11 @@ gfortran -c -O3 -fopenmp prtarc.f gfortran -c -O3 -fopenmp prtcif.f gfortran -c -O3 -fopenmp prtdyn.f gfortran -c -O3 -fopenmp prterr.f -gfortran -c -O3 -fopenmp prtfrc.f gfortran -c -O3 -fopenmp prtint.f gfortran -c -O3 -fopenmp prtmol2.f gfortran -c -O3 -fopenmp prtpdb.f gfortran -c -O3 -fopenmp prtprm.f gfortran -c -O3 -fopenmp prtseq.f -gfortran -c -O3 -fopenmp prtuind.f -gfortran -c -O3 -fopenmp prtvel.f gfortran -c -O3 -fopenmp prtxyz.f gfortran -c -O3 -fopenmp pss.f gfortran -c -O3 -fopenmp pssrigid.f diff --git a/macos/gfortran/debug.make b/macos/gfortran/debug.make index e62df213a..070640104 100755 --- a/macos/gfortran/debug.make +++ b/macos/gfortran/debug.make @@ -528,14 +528,11 @@ gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prt gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtcif.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtdyn.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prterr.f -gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtfrc.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtint.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtmol2.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtpdb.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtprm.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtseq.f -gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtuind.f -gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtvel.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized prtxyz.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized pss.f gfortran -c -Og -g -fbacktrace -fcheck=bounds -Wunused -Wmaybe-uninitialized pssrigid.f diff --git a/macos/gfortran/library.make b/macos/gfortran/library.make index b6846b1ee..3a6152b11 100755 --- a/macos/gfortran/library.make +++ b/macos/gfortran/library.make @@ -439,14 +439,11 @@ prtarc.o \ prtcif.o \ prtdyn.o \ prterr.o \ -prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ -prtuind.o \ -prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ diff --git a/make/Makefile b/make/Makefile index 87b3ba33a..f370d48ab 100644 --- a/make/Makefile +++ b/make/Makefile @@ -592,14 +592,11 @@ OBJS = action.o \ prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ pss.o \ pssrigid.o \ @@ -1524,14 +1521,11 @@ libtinker.a: ${OBJS} prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ @@ -2106,14 +2100,11 @@ protein.o: atomid.o atoms.o couple.o files.o group.o inform.o iounit.o katoms.o prtarc.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o output.o titles.o usage.o prtcif.o: ascii.o bound.o boxes.o files.o pdb.o sequen.o titles.o prtdyn.o: atoms.o boxes.o files.o group.o mdstuf.o moldyn.o rgddyn.o titles.o -prtfrc.o: atomid.o atoms.o bound.o boxes.o couple.o deriv.o files.o inform.o titles.o prtint.o: atomid.o atoms.o files.o inform.o titles.o zclose.o zcoord.o prtmol2.o: angbnd.o atmlst.o atomid.o atoms.o bndstr.o couple.o files.o iounit.o ptable.o ring.o titles.o tors.o prtpdb.o: bound.o boxes.o files.o pdb.o sequen.o titles.o prtprm.o: angpot.o bndpot.o chgpot.o fields.o kanang.o kangs.o kantor.o katoms.o kbonds.o kcflux.o kchrge.o kcpen.o kctrn.o kdipol.o kdsp.o kexpl.o khbond.o kiprop.o kitors.o kmulti.o kopbnd.o kopdst.o korbs.o kpitor.o kpolpr.o kpolr.o krepl.o ksolut.o kstbnd.o ksttor.o ktorsn.o ktrtor.o kurybr.o kvdwpr.o kvdws.o mplpot.o polpot.o sizes.o urypot.o vdwpot.o prtseq.o: files.o sequen.o -prtuind.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o polar.o potent.o solpot.o titles.o units.o -prtvel.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o moldyn.o titles.o prtxyz.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o titles.o pss.o: atoms.o files.o hescut.o inform.o iounit.o math.o omega.o refer.o tree.o warp.o zcoord.o pssrigid.o: atoms.o files.o group.o inform.o iounit.o math.o minima.o molcul.o refer.o rigid.o sizes.o warp.o diff --git a/make/Makefile-apbs b/make/Makefile-apbs index 9ffe1d84d..1f62207cd 100644 --- a/make/Makefile-apbs +++ b/make/Makefile-apbs @@ -613,14 +613,11 @@ OBJS = action.o \ prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ pss.o \ pssrigid.o \ @@ -1548,14 +1545,11 @@ libtinker.a: ${OBJS} prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ @@ -2130,14 +2124,11 @@ protein.o: atomid.o atoms.o couple.o files.o group.o inform.o iounit.o katoms.o prtarc.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o output.o titles.o usage.o prtcif.o: ascii.o bound.o boxes.o files.o pdb.o sequen.o titles.o prtdyn.o: atoms.o boxes.o files.o group.o mdstuf.o moldyn.o rgddyn.o titles.o -prtfrc.o: atomid.o atoms.o bound.o boxes.o couple.o deriv.o files.o inform.o titles.o prtint.o: atomid.o atoms.o files.o inform.o titles.o zclose.o zcoord.o prtmol2.o: angbnd.o atmlst.o atomid.o atoms.o bndstr.o couple.o files.o iounit.o ptable.o ring.o titles.o tors.o prtpdb.o: bound.o boxes.o files.o pdb.o sequen.o titles.o prtprm.o: angpot.o bndpot.o chgpot.o fields.o kanang.o kangs.o kantor.o katoms.o kbonds.o kcflux.o kchrge.o kcpen.o kctrn.o kdipol.o kdsp.o kexpl.o khbond.o kiprop.o kitors.o kmulti.o kopbnd.o kopdst.o korbs.o kpitor.o kpolpr.o kpolr.o krepl.o ksolut.o kstbnd.o ksttor.o ktorsn.o ktrtor.o kurybr.o kvdwpr.o kvdws.o mplpot.o polpot.o sizes.o urypot.o vdwpot.o prtseq.o: files.o sequen.o -prtuind.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o polar.o potent.o solpot.o titles.o units.o -prtvel.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o moldyn.o titles.o prtxyz.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o titles.o pss.o: atoms.o files.o hescut.o inform.o iounit.o math.o omega.o refer.o tree.o warp.o zcoord.o pssrigid.o: atoms.o files.o group.o inform.o iounit.o math.o minima.o molcul.o refer.o rigid.o sizes.o warp.o diff --git a/make/Makefile-ffe-linux b/make/Makefile-ffe-linux index 6eaeca4ec..2674efc1f 100644 --- a/make/Makefile-ffe-linux +++ b/make/Makefile-ffe-linux @@ -640,14 +640,11 @@ OBJS = action.o \ prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ pss.o \ pssrigid.o \ @@ -1575,14 +1572,11 @@ libtinker.a: ${OBJS} prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ @@ -2157,14 +2151,11 @@ protein.o: atomid.o atoms.o couple.o files.o group.o inform.o iounit.o katoms.o prtarc.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o output.o titles.o usage.o prtcif.o: ascii.o bound.o boxes.o files.o pdb.o sequen.o titles.o prtdyn.o: atoms.o boxes.o files.o group.o mdstuf.o moldyn.o rgddyn.o titles.o -prtfrc.o: atomid.o atoms.o bound.o boxes.o couple.o deriv.o files.o inform.o titles.o prtint.o: atomid.o atoms.o files.o inform.o titles.o zclose.o zcoord.o prtmol2.o: angbnd.o atmlst.o atomid.o atoms.o bndstr.o couple.o files.o iounit.o ptable.o ring.o titles.o tors.o prtpdb.o: bound.o boxes.o files.o pdb.o sequen.o titles.o prtprm.o: angpot.o bndpot.o chgpot.o fields.o kanang.o kangs.o kantor.o katoms.o kbonds.o kcflux.o kchrge.o kcpen.o kctrn.o kdipol.o kdsp.o kexpl.o khbond.o kiprop.o kitors.o kmulti.o kopbnd.o kopdst.o korbs.o kpitor.o kpolpr.o kpolr.o krepl.o ksolut.o kstbnd.o ksttor.o ktorsn.o ktrtor.o kurybr.o kvdwpr.o kvdws.o mplpot.o polpot.o sizes.o urypot.o vdwpot.o prtseq.o: files.o sequen.o -prtuind.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o polar.o potent.o solpot.o titles.o units.o -prtvel.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o moldyn.o titles.o prtxyz.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o titles.o pss.o: atoms.o files.o hescut.o inform.o iounit.o math.o omega.o refer.o tree.o warp.o zcoord.o pssrigid.o: atoms.o files.o group.o inform.o iounit.o math.o minima.o molcul.o refer.o rigid.o sizes.o warp.o diff --git a/make/Makefile-ffe-macos b/make/Makefile-ffe-macos index d843cab45..664fd90ad 100644 --- a/make/Makefile-ffe-macos +++ b/make/Makefile-ffe-macos @@ -640,14 +640,11 @@ OBJS = action.o \ prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ pss.o \ pssrigid.o \ @@ -1575,14 +1572,11 @@ libtinker.a: ${OBJS} prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ @@ -2157,14 +2151,11 @@ protein.o: atomid.o atoms.o couple.o files.o group.o inform.o iounit.o katoms.o prtarc.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o output.o titles.o usage.o prtcif.o: ascii.o bound.o boxes.o files.o pdb.o sequen.o titles.o prtdyn.o: atoms.o boxes.o files.o group.o mdstuf.o moldyn.o rgddyn.o titles.o -prtfrc.o: atomid.o atoms.o bound.o boxes.o couple.o deriv.o files.o inform.o titles.o prtint.o: atomid.o atoms.o files.o inform.o titles.o zclose.o zcoord.o prtmol2.o: angbnd.o atmlst.o atomid.o atoms.o bndstr.o couple.o files.o iounit.o ptable.o ring.o titles.o tors.o prtpdb.o: bound.o boxes.o files.o pdb.o sequen.o titles.o prtprm.o: angpot.o bndpot.o chgpot.o fields.o kanang.o kangs.o kantor.o katoms.o kbonds.o kcflux.o kchrge.o kcpen.o kctrn.o kdipol.o kdsp.o kexpl.o khbond.o kiprop.o kitors.o kmulti.o kopbnd.o kopdst.o korbs.o kpitor.o kpolpr.o kpolr.o krepl.o ksolut.o kstbnd.o ksttor.o ktorsn.o ktrtor.o kurybr.o kvdwpr.o kvdws.o mplpot.o polpot.o sizes.o urypot.o vdwpot.o prtseq.o: files.o sequen.o -prtuind.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o polar.o potent.o solpot.o titles.o units.o -prtvel.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o moldyn.o titles.o prtxyz.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o titles.o pss.o: atoms.o files.o hescut.o inform.o iounit.o math.o omega.o refer.o tree.o warp.o zcoord.o pssrigid.o: atoms.o files.o group.o inform.o iounit.o math.o minima.o molcul.o refer.o rigid.o sizes.o warp.o diff --git a/make/Makefile-ffe-windows b/make/Makefile-ffe-windows index 43e2415a5..11e06d03c 100644 --- a/make/Makefile-ffe-windows +++ b/make/Makefile-ffe-windows @@ -640,14 +640,11 @@ OBJS = action.o \ prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ pss.o \ pssrigid.o \ @@ -1584,14 +1581,11 @@ libtinker.a: ${OBJS} prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ @@ -2166,14 +2160,11 @@ protein.o: atomid.o atoms.o couple.o files.o group.o inform.o iounit.o katoms.o prtarc.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o output.o titles.o usage.o prtcif.o: ascii.o bound.o boxes.o files.o pdb.o sequen.o titles.o prtdyn.o: atoms.o boxes.o files.o group.o mdstuf.o moldyn.o rgddyn.o titles.o -prtfrc.o: atomid.o atoms.o bound.o boxes.o couple.o deriv.o files.o inform.o titles.o prtint.o: atomid.o atoms.o files.o inform.o titles.o zclose.o zcoord.o prtmol2.o: angbnd.o atmlst.o atomid.o atoms.o bndstr.o couple.o files.o iounit.o ptable.o ring.o titles.o tors.o prtpdb.o: bound.o boxes.o files.o pdb.o sequen.o titles.o prtprm.o: angpot.o bndpot.o chgpot.o fields.o kanang.o kangs.o kantor.o katoms.o kbonds.o kcflux.o kchrge.o kcpen.o kctrn.o kdipol.o kdsp.o kexpl.o khbond.o kiprop.o kitors.o kmulti.o kopbnd.o kopdst.o korbs.o kpitor.o kpolpr.o kpolr.o krepl.o ksolut.o kstbnd.o ksttor.o ktorsn.o ktrtor.o kurybr.o kvdwpr.o kvdws.o mplpot.o polpot.o sizes.o urypot.o vdwpot.o prtseq.o: files.o sequen.o -prtuind.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o polar.o potent.o solpot.o titles.o units.o -prtvel.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o moldyn.o titles.o prtxyz.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o titles.o pss.o: atoms.o files.o hescut.o inform.o iounit.o math.o omega.o refer.o tree.o warp.o zcoord.o pssrigid.o: atoms.o files.o group.o inform.o iounit.o math.o minima.o molcul.o refer.o rigid.o sizes.o warp.o diff --git a/make/depend.make b/make/depend.make index b095113ed..c08e79875 100755 --- a/make/depend.make +++ b/make/depend.make @@ -471,14 +471,11 @@ document 5 protein.f | grep 'o:' document 5 prtarc.f | grep 'o:' document 5 prtcif.f | grep 'o:' document 5 prtdyn.f | grep 'o:' -document 5 prtfrc.f | grep 'o:' document 5 prtint.f | grep 'o:' document 5 prtmol2.f | grep 'o:' document 5 prtpdb.f | grep 'o:' document 5 prtprm.f | grep 'o:' document 5 prtseq.f | grep 'o:' -document 5 prtuind.f | grep 'o:' -document 5 prtvel.f | grep 'o:' document 5 prtxyz.f | grep 'o:' document 5 pss.f | grep 'o:' document 5 pssrigid.f | grep 'o:' diff --git a/openmm/Makefile b/openmm/Makefile index 16600fb3a..cf40962f2 100644 --- a/openmm/Makefile +++ b/openmm/Makefile @@ -622,14 +622,11 @@ OBJS = action.o \ prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ pss.o \ pssrigid.o \ @@ -1256,14 +1253,11 @@ libtinker.a: ${OBJS} prtcif.o \ prtdyn.o \ prterr.o \ - prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ - prtuind.o \ - prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ @@ -1845,14 +1839,11 @@ protein.o: atomid.o atoms.o couple.o files.o group.o inform.o iounit.o katoms.o prtarc.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o output.o titles.o usage.o prtcif.o: ascii.o bound.o boxes.o files.o pdb.o sequen.o titles.o prtdyn.o: atoms.o boxes.o files.o group.o mdstuf.o moldyn.o rgddyn.o titles.o -prtfrc.o: atomid.o atoms.o bound.o boxes.o couple.o deriv.o files.o inform.o titles.o prtint.o: atomid.o atoms.o files.o inform.o titles.o zclose.o zcoord.o prtmol2.o: angbnd.o atmlst.o atomid.o atoms.o bndstr.o couple.o files.o iounit.o ptable.o ring.o titles.o tors.o prtpdb.o: bound.o boxes.o files.o pdb.o sequen.o titles.o prtprm.o: angpot.o bndpot.o chgpot.o fields.o kanang.o kangs.o kantor.o katoms.o kbonds.o kcflux.o kchrge.o kcpen.o kctrn.o kdipol.o kdsp.o kexpl.o khbond.o kiprop.o kitors.o kmulti.o kopbnd.o kopdst.o korbs.o kpitor.o kpolpr.o kpolr.o krepl.o ksolut.o kstbnd.o ksttor.o ktorsn.o ktrtor.o kurybr.o kvdwpr.o kvdws.o mplpot.o polpot.o sizes.o urypot.o vdwpot.o prtseq.o: files.o sequen.o -prtuind.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o polar.o potent.o solpot.o titles.o units.o -prtvel.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o moldyn.o titles.o prtxyz.o: atomid.o atoms.o bound.o boxes.o couple.o files.o inform.o titles.o pss.o: atoms.o files.o hescut.o inform.o iounit.o math.o omega.o refer.o tree.o warp.o zcoord.o pssrigid.o: atoms.o files.o group.o inform.o iounit.o math.o minima.o molcul.o refer.o rigid.o sizes.o warp.o diff --git a/source/active.f b/source/active.f index f8bfdae2e..c9244695d 100644 --- a/source/active.f +++ b/source/active.f @@ -243,3 +243,228 @@ subroutine active deallocate (fixed) return end +c +c +c ################################################################# +c ## ## +c ## subroutine saveonly -- set the list of save coord atoms ## +c ## ## +c ################################################################# +c +c +c "saveonly" sets the list of atoms that are used during +c coordinate saving routines +c +c + subroutine saveonly + use atoms + use iounit + use keys + use output + implicit none + integer i,j,next + integer nfixed + integer, allocatable :: fixed(:) + character*20 keyword + character*240 record + character*240 string + logical header + logical, allocatable :: saved(:) +c +c +c perform dynamic allocation of some global arrays +c + if (allocated(ionly)) deallocate (ionly) + if (allocated(ionlyinv)) deallocate (ionlyinv) + if (allocated(print3n)) deallocate (print3n) + allocate (ionly(n)) + allocate (ionlyinv(n)) + allocate (print3n(3,n)) +c +c perform dynamic allocation of some local arrays +c + allocate (fixed(n)) + allocate (saved(n)) +c +c set defaults for the numbers and lists of saved atoms +c + onlysave = .false. + nonly = 0 + do i = 1, n + ionly(i) = 0 + ionlyinv(i) = 0 + end do + nfixed = 0 + do i = 1, n + fixed(i) = 0 + saved(i) = .false. + end do +c +c get any keywords containing save-only atom parameters +c + do j = 1, nkey + next = 1 + record = keyline(j) + call gettext (record,keyword,next) + call upcase (keyword) + string = record(next:240) +c +c get any lists of atoms whose coordinates should be saved +c + if (keyword(1:10) .eq. 'SAVE-ONLY ') then + read (string,*,err=10,end=10) (fixed(i),i=nfixed+1,n) + 10 continue + do while (fixed(nfixed+1) .ne. 0) + nfixed = nfixed + 1 + end do + end if + end do +c +c remove saved atoms not in the system +c + header = .true. + do i = 1, n + if (abs(fixed(i)) .gt. n) then + fixed(i) = 0 + if (header) then + header = .false. + write (iout,20) + 20 format (/,' SAVEONLY -- Warning, Illegal Atom Number', + & ' in SAVE-ONLY Atom List') + end if + end if + end do +c +c set saved atoms to only those marked as save +c + i = 1 + do while (fixed(i) .ne. 0) + if (fixed(i) .gt. 0) then + j = fixed(i) + saved(j) = .true. + i = i + 1 + else + do j = abs(fixed(i)), abs(fixed(i+1)) + saved(j) = .true. + end do + i = i + 2 + end if + end do + do i = 1, n + if (saved(i)) then + nonly = nonly + 1 + ionly(nonly) = i + ionlyinv(i) = nonly + end if + end do + if (nonly > 0) onlysave = .true. +c +c perform deallocation of some local arrays +c + deallocate (fixed) + deallocate (saved) + return + end +c +c +c ################################################################## +c ## ## +c ## subroutine msystem -- set exclusion for moment of system ## +c ## ## +c ################################################################## +c +c +c "msystem" sets the list of atoms that are excluded while +c computing moment of system +c +c + subroutine msystem + use atoms + use iounit + use keys + use moment + implicit none + integer i,j,next + integer nfixed + integer, allocatable :: fixed(:) + character*20 keyword + character*240 record + character*240 string + logical header +c +c +c perform dynamic allocation of some global arrays +c + if (allocated(momuse)) deallocate (momuse) + allocate (momuse(n)) +c +c perform dynamic allocation of some local arrays +c + allocate (fixed(n)) +c +c set defaults for the numbers and lists of atoms to be used +c + do i = 1, n + momuse(i) = .true. + end do + nfixed = 0 + do i = 1, n + fixed(i) = 0 + end do +c +c get any keywords containing exc-moment atom parameters +c + do j = 1, nkey + next = 1 + record = keyline(j) + call gettext (record,keyword,next) + call upcase (keyword) + string = record(next:240) +c +c get any lists of atoms whose coordinates should be used +c + if (keyword(1:13) .eq. 'EXC-MOMENT ') then + read (string,*,err=10,end=10) (fixed(i),i=nfixed+1,n) + 10 continue + do while (fixed(nfixed+1) .ne. 0) + nfixed = nfixed + 1 + end do + end if + end do +c +c remove used atoms not in the system +c + header = .true. + do i = 1, n + if (abs(fixed(i)) .gt. n) then + fixed(i) = 0 + if (header) then + header = .false. + write (iout,20) + 20 format (/,' MSYSTEM -- Warning, Illegal Atom Number', + & ' in EXC-MOMENT Atom List') + end if + end if + end do +c +c set inactive atoms to false +c + i = 1 + do while (fixed(i) .ne. 0) + if (fixed(i) .gt. 0) then + j = fixed(i) + momuse(j) = .false. + i = i + 1 + else + do j = abs(fixed(i)), abs(fixed(i+1)) + momuse(j) = .false. + end do + i = i + 2 + end if + end do +c +c perform deallocation of some local arrays +c + deallocate (fixed) + return + end diff --git a/source/boxes.f b/source/boxes.f index 9eb1e38b4..18d8eabaf 100644 --- a/source/boxes.f +++ b/source/boxes.f @@ -31,6 +31,9 @@ c gamma_cos cosine of the gamma periodic box angle c beta_term term used in generating triclinic box c gamma_term term used in generating triclinic box +c xcenter x-coordinate of center of mass of system in Angstroms +c ycenter y-coordinate of center of mass of system in Angstroms +c zcenter z-coordinate of center of mass of system in Angstroms c lvec real space lattice vectors as matrix rows c recip reciprocal lattice vectors as matrix columns c orthogonal flag to mark periodic box as orthogonal @@ -57,6 +60,7 @@ module boxes real*8 gamma_cos real*8 beta_term real*8 gamma_term + real*8 xcenter,ycenter,zcenter real*8 lvec(3,3) real*8 recip(3,3) logical orthogonal diff --git a/source/center.f b/source/center.f index 48c29fd36..0c9b08797 100644 --- a/source/center.f +++ b/source/center.f @@ -76,3 +76,58 @@ subroutine center (n1,x1,y1,z1,n2,x2,y2,z2,xmid,ymid,zmid) end do return end +c +c +c ####################################################### +c ## ## +c ## subroutine compcent -- compute center of mass ## +c ## ## +c ####################################################### +c +c +c "compcent" computes the center of mass +c +c + subroutine compcent (xcm,ycm,zcm) + use atomid + use atoms + use boxes + implicit none + integer i + real*8 weigh + real*8 xcm,ycm,zcm +c +c +c find the center of mass of the set of active atoms +c + weigh = 0.0d0 + xcenter = 0.0d0 + ycenter = 0.0d0 + zcenter = 0.0d0 +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(n,x,y,z,xcenter,ycenter,zcenter,weigh,mass) +!$OMP DO reduction(+:xcenter,ycenter,zcenter,weigh) + do i = 1, n + weigh = weigh + mass(i) + xcenter = xcenter + x(i)*mass(i) + ycenter = ycenter + y(i)*mass(i) + zcenter = zcenter + z(i)*mass(i) + end do +!$OMP END DO +!$OMP END PARALLEL + if (weigh .ne. 0.0d0) then + xcenter = xcenter / weigh + ycenter = ycenter / weigh + zcenter = zcenter / weigh + end if +c +c copy xcenter, ycenter, zcenter to xcm, ycm, zcm +c + xcm = xcenter + ycm = ycenter + zcm = zcenter + return + end diff --git a/source/control.f b/source/control.f index 4458f9177..25778cae5 100644 --- a/source/control.f +++ b/source/control.f @@ -41,6 +41,8 @@ subroutine control cyclesave = .false. noversion = .false. overwrite = .false. + coordsave = .true. + dynsave = .true. c c check for control parameters on the command line c @@ -90,6 +92,10 @@ subroutine control noversion = .true. else if (keyword(1:10) .eq. 'OVERWRITE ') then overwrite = .true. + else if (keyword(1:8) .eq. 'NOCOORD ') then + coordsave = .false. + else if (keyword(1:6) .eq. 'NODYN ') then + dynsave = .false. end if 10 continue end do diff --git a/source/dynamic.f b/source/dynamic.f index 5708186cc..5d343501d 100644 --- a/source/dynamic.f +++ b/source/dynamic.f @@ -22,17 +22,20 @@ program dynamic use bath use bndstr use bound + use extfld use inform use iounit use keys use mdstuf use potent use stodyn + use output use usage implicit none integer i,next,mode integer istep,nstep real*8 dt,dtsave + real*8 phs logical exist character*20 keyword character*240 record @@ -226,55 +229,67 @@ program dynamic c call mdinit (dt) c +c only allow Montecarlo or anisotropic barostat +c for NPT + extfield simulation +c + if (use_exfld .and. mode.eq.4) then + if (barostat.ne.'MONTECARLO') then + write (iout,340) + 340 format (/,' DYNAMIC -- NPT with External Field Should', + & ' Use MonteCarlo Barostat') + call fatal + end if + end if +c c print out a header line for the dynamics computation c if (integrate .eq. 'VERLET') then - write (iout,340) - 340 format (/,' Molecular Dynamics Trajectory via', - & ' Velocity Verlet Algorithm') - else if (integrate .eq. 'BEEMAN') then write (iout,350) 350 format (/,' Molecular Dynamics Trajectory via', + & ' Velocity Verlet Algorithm') + else if (integrate .eq. 'BEEMAN') then + write (iout,360) + 360 format (/,' Molecular Dynamics Trajectory via', & ' Modified Beeman Algorithm') else if (integrate .eq. 'BAOAB') then - write (iout,360) - 360 format (/,' Stochastic Dynamics Trajectory via', - & ' BAOAB Algorithm') - else if (integrate .eq. 'OBABO') then write (iout,370) 370 format (/,' Stochastic Dynamics Trajectory via', + & ' BAOAB Algorithm') + else if (integrate .eq. 'OBABO') then + write (iout,380) + 380 format (/,' Stochastic Dynamics Trajectory via', & ' OBABO Algorithm') else if (integrate .eq. 'NOSE-HOOVER') then - write (iout,380) - 380 format (/,' Molecular Dynamics Trajectory via', + write (iout,390) + 390 format (/,' Molecular Dynamics Trajectory via', & ' Nose-Hoover NPT Algorithm') else if (integrate .eq. 'STOCHASTIC') then - write (iout,390) - 390 format (/,' Stochastic Dynamics Trajectory via', - & ' Velocity Verlet Algorithm') - else if (integrate .eq. 'GHMC') then write (iout,400) 400 format (/,' Stochastic Dynamics Trajectory via', + & ' Velocity Verlet Algorithm') + else if (integrate .eq. 'GHMC') then + write (iout,410) + 410 format (/,' Stochastic Dynamics Trajectory via', & ' Generalized Hybrid Monte Carlo') else if (integrate .eq. 'RIGIDBODY') then - write (iout,410) - 410 format (/,' Molecular Dynamics Trajectory via', - & ' Rigid Body Algorithm') - else if (integrate .eq. 'VRESPA') then write (iout,420) 420 format (/,' Molecular Dynamics Trajectory via', - & ' Verlet r-RESPA MTS Algorithm') - else if (integrate .eq. 'BRESPA') then + & ' Rigid Body Algorithm') + else if (integrate .eq. 'VRESPA') then write (iout,430) 430 format (/,' Molecular Dynamics Trajectory via', - & ' Beeman r-RESPA MTS Algorithm') - else if (integrate .eq. 'SRESPA') then + & ' Verlet r-RESPA MTS Algorithm') + else if (integrate .eq. 'BRESPA') then write (iout,440) 440 format (/,' Molecular Dynamics Trajectory via', - & ' BAOAB r-RESPA MTS Algorithm') - else + & ' Beeman r-RESPA MTS Algorithm') + else if (integrate .eq. 'SRESPA') then write (iout,450) 450 format (/,' Molecular Dynamics Trajectory via', + & ' BAOAB r-RESPA MTS Algorithm') + else + write (iout,460) + 460 format (/,' Molecular Dynamics Trajectory via', & ' Modified Beeman Algorithm') end if flush (iout) @@ -282,6 +297,12 @@ program dynamic c integrate equations of motion to take a time step c do istep = 1, nstep + if (use_exfld .and. use_exfreq) then + phs = sin(exfreq * dble(istep-1) * dt) + do i = 1, 3 + texfld(i) = phs * exfld(i) + end do + end if if (integrate .eq. 'VERLET') then call verlet (istep,dt) else if (integrate .eq. 'BEEMAN') then @@ -309,6 +330,10 @@ program dynamic end if end do c +c save dynamic at the end if it was not saved during simulation +c + if (.not. dynsave) call prtdyn +c c perform any final tasks before program exit c call final diff --git a/source/epolar.f b/source/epolar.f index b85ace989..69b0f748c 100644 --- a/source/epolar.f +++ b/source/epolar.f @@ -483,7 +483,7 @@ subroutine epolar0a do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e end do @@ -594,7 +594,7 @@ subroutine epolar0b !$OMP& n13,i13,n14,i14,n15,i15,np11,ip11,np12,ip12,np13,ip13,np14,ip14, !$OMP& p2scale,p3scale,p4scale,p5scale,p2iscale,p3iscale,p4iscale, !$OMP& p5iscale,nelst,elst,use_thole,use_chgpen,use_bounds,f,off2, -!$OMP& exfld,use_exfld) +!$OMP& texfld,use_exfld) !$OMP& firstprivate(pscale) shared (ep) !$OMP DO reduction(+:ep) c @@ -774,7 +774,7 @@ subroutine epolar0b do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e end do @@ -1360,7 +1360,7 @@ subroutine epreal0c do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e end do @@ -1568,7 +1568,7 @@ subroutine epreal0d !$OMP& n13,i13,n14,i14,n15,i15,np11,ip11,np12,ip12,np13,ip13,np14,ip14, !$OMP& p2scale,p3scale,p4scale,p5scale,p2iscale,p3iscale,p4iscale, !$OMP& p5iscale,nelst,elst,use_thole,use_chgpen,use_bounds,off2,f, -!$OMP& exfld,use_exfld) +!$OMP& texfld,use_exfld) !$OMP& firstprivate(pscale) shared (ep) !$OMP DO reduction(+:ep) c @@ -1765,7 +1765,7 @@ subroutine epreal0d do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e end do diff --git a/source/epolar3.f b/source/epolar3.f index f2cbf2268..4490810f2 100644 --- a/source/epolar3.f +++ b/source/epolar3.f @@ -556,7 +556,7 @@ subroutine epolar3a do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e nep = nep + 1 @@ -692,7 +692,7 @@ subroutine epolar3b !$OMP& n13,i13,n14,i14,n15,i15,np11,ip11,np12,ip12,np13,ip13,np14,ip14, !$OMP& p2scale,p3scale,p4scale,p5scale,p2iscale,p3iscale,p4iscale, !$OMP& p5iscale,nelst,elst,use_thole,use_chgpen,use_bounds,off2,f, -!$OMP& exfld,use_exfld,molcule,name,verbose,debug,header,iout) +!$OMP& texfld,molcule,name,verbose,debug,header,iout,use_exfld) !$OMP& firstprivate(pscale) shared (ep,nep,aep,einter) !$OMP DO reduction(+:ep,nep,aep,einter) c @@ -897,7 +897,7 @@ subroutine epolar3b do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e nep = nep + 1 @@ -1588,7 +1588,7 @@ subroutine epreal3c do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e nep = nep + 1 @@ -1832,7 +1832,7 @@ subroutine epreal3d !$OMP& n13,i13,n14,i14,n15,i15,np11,ip11,np12,ip12,np13,ip13,np14,ip14, !$OMP& p2scale,p3scale,p4scale,p5scale,p2iscale,p3iscale,p4iscale, !$OMP& p5iscale,nelst,elst,use_thole,use_chgpen,use_bounds,off2,f, -!$OMP& exfld,use_exfld,molcule,name,verbose,debug,header,iout) +!$OMP& texfld,molcule,name,verbose,debug,header,iout,use_exfld) !$OMP& firstprivate(pscale) shared (ep,nep,aep,einter) !$OMP DO reduction(+:ep,nep,aep,einter) c @@ -2062,7 +2062,7 @@ subroutine epreal3d do i = 1, npole e = 0.0d0 do j = 1, 3 - e = e - f*uind(j,i)*exfld(j) + e = e - f*uind(j,i)*texfld(j) end do ep = ep + e nep = nep + 1 diff --git a/source/exfield.f b/source/exfield.f index 8266ed667..4ad792d47 100644 --- a/source/exfield.f +++ b/source/exfield.f @@ -18,6 +18,7 @@ c subroutine exfield (mode,exf) use atoms + use bound use charge use chgpot use energi @@ -37,11 +38,15 @@ subroutine exfield (mode,exf) exf = 0.0d0 f = electric / dielec c +c maintain any periodic boundary conditions +c + if (use_bounds) call bounds +c c calculate external field energy over partial charges c if (mode .eq. 'CHARGE') then !$OMP PARALLEL default(private) shared(nion,iion,use, -!$OMP& x,y,z,f,pchg,exfld,exf) +!$OMP& x,y,z,f,pchg,texfld,exf) !$OMP DO reduction(+:exf) do ii = 1, nion i = iion(ii) @@ -50,7 +55,7 @@ subroutine exfield (mode,exf) yi = y(i) zi = z(i) ci = pchg(i) - phi = xi*exfld(1) + yi*exfld(2) + zi*exfld(3) + phi = xi*texfld(1) + yi*texfld(2) + zi*texfld(3) e = -f * ci * phi exf = exf + e end if @@ -63,7 +68,7 @@ subroutine exfield (mode,exf) c if (mode .eq. 'MPOLE') then !$OMP PARALLEL default(private) shared(npole,ipole,use, -!$OMP& x,y,z,f,rpole,exfld,exf) +!$OMP& x,y,z,f,rpole,texfld,exf) !$OMP DO reduction(+:exf) do ii = 1, npole i = ipole(ii) @@ -72,12 +77,12 @@ subroutine exfield (mode,exf) yi = y(i) zi = z(i) ci = rpole(1,i) - phi = xi*exfld(1) + yi*exfld(2) + zi*exfld(3) + phi = xi*texfld(1) + yi*texfld(2) + zi*texfld(3) dix = rpole(2,i) diy = rpole(3,i) diz = rpole(4,i) - e = -f * (ci*phi + dix*exfld(1) - & + diy*exfld(2) + diz*exfld(3)) + e = -f * (ci*phi + dix*texfld(1) + & + diy*texfld(2) + diz*texfld(3)) exf = exf + e end if end do @@ -101,6 +106,7 @@ subroutine exfield (mode,exf) c subroutine exfield1 (mode,exf) use atoms + use bound use charge use chgpot use deriv @@ -131,11 +137,15 @@ subroutine exfield1 (mode,exf) exf = 0.0d0 f = electric / dielec c +c maintain any periodic boundary conditions +c + if (use_bounds) call bounds +c c calculate energy and derivatives over partial charges c if (mode .eq. 'CHARGE') then !$OMP PARALLEL default(private) shared(nion,iion,use, -!$OMP& x,y,z,f,pchg,exfld,exf,dec,vir) +!$OMP& x,y,z,f,pchg,texfld,exf,dec,vir) !$OMP DO reduction(+:exf,dec,vir) do ii = 1, nion i = iion(ii) @@ -144,15 +154,15 @@ subroutine exfield1 (mode,exf) yi = y(i) zi = z(i) ci = pchg(i) - phi = xi*exfld(1) + yi*exfld(2) + zi*exfld(3) + phi = xi*texfld(1) + yi*texfld(2) + zi*texfld(3) e = -f * ci * phi exf = exf + e c c gradient and virial components from charge interactions c - frx = -f * exfld(1) * ci - fry = -f * exfld(2) * ci - frz = -f * exfld(3) * ci + frx = -f * texfld(1) * ci + fry = -f * texfld(2) * ci + frz = -f * texfld(3) * ci dec(1,i) = dec(1,i) + frx dec(2,i) = dec(2,i) + fry dec(3,i) = dec(3,i) + frz @@ -184,7 +194,7 @@ subroutine exfield1 (mode,exf) c if (mode .eq. 'MPOLE') then !$OMP PARALLEL default(private) shared(npole,ipole,use, -!$OMP& x,y,z,xaxis,yaxis,zaxis,f,rpole,exfld,exf,dem,vir) +!$OMP& x,y,z,xaxis,yaxis,zaxis,f,rpole,texfld,exf,dem,vir) !$OMP DO reduction(+:exf,dem,vir) do ii = 1, npole i = ipole(ii) @@ -196,16 +206,16 @@ subroutine exfield1 (mode,exf) dix = rpole(2,i) diy = rpole(3,i) diz = rpole(4,i) - phi = xi*exfld(1) + yi*exfld(2) + zi*exfld(3) - e = -f * (ci*phi + dix*exfld(1) - & + diy*exfld(2) + diz*exfld(3)) + phi = xi*texfld(1) + yi*texfld(2) + zi*texfld(3) + e = -f * (ci*phi + dix*texfld(1) + & + diy*texfld(2) + diz*texfld(3)) exf = exf + e c c gradient and virial components from dipole interactions c - tem(1) = f * (diy*exfld(3)-diz*exfld(2)) - tem(2) = f * (diz*exfld(1)-dix*exfld(3)) - tem(3) = f * (dix*exfld(2)-diy*exfld(1)) + tem(1) = f * (diy*texfld(3)-diz*texfld(2)) + tem(2) = f * (diz*texfld(1)-dix*texfld(3)) + tem(3) = f * (dix*texfld(2)-diy*texfld(1)) call torque (i,tem,fix,fiy,fiz,dem) iz = zaxis(i) ix = xaxis(i) @@ -234,9 +244,9 @@ subroutine exfield1 (mode,exf) c c gradient and virial components from monopole interactions c - frx = -f * exfld(1) * ci - fry = -f * exfld(2) * ci - frz = -f * exfld(3) * ci + frx = -f * texfld(1) * ci + fry = -f * texfld(2) * ci + frz = -f * texfld(3) * ci dem(1,i) = dem(1,i) + frx dem(2,i) = dem(2,i) + fry dem(3,i) = dem(3,i) + frz @@ -283,6 +293,7 @@ subroutine exfield3 (mode,exf) use action use analyz use atoms + use bound use charge use chgpot use energi @@ -302,11 +313,15 @@ subroutine exfield3 (mode,exf) exf = 0.0d0 f = electric / dielec c +c maintain any periodic boundary conditions +c + if (use_bounds) call bounds +c c calculate energy and partitioning over partial charges c if (mode .eq. 'CHARGE') then !$OMP PARALLEL default(private) shared(nion,iion,use, -!$OMP& x,y,z,f,pchg,exfld,exf,nec,aec) +!$OMP& x,y,z,f,pchg,texfld,exf,nec,aec) !$OMP DO reduction(+:exf,nec,aec) do ii = 1, nion i = iion(ii) @@ -315,7 +330,7 @@ subroutine exfield3 (mode,exf) yi = y(i) zi = z(i) ci = pchg(i) - phi = xi*exfld(1) + yi*exfld(2) + zi*exfld(3) + phi = xi*texfld(1) + yi*texfld(2) + zi*texfld(3) e = -f * ci * phi exf = exf + e nec = nec + 1 @@ -330,7 +345,7 @@ subroutine exfield3 (mode,exf) c if (mode .eq. 'MPOLE') then !$OMP PARALLEL default(private) shared(npole,ipole,use, -!$OMP& x,y,z,f,rpole,exfld,exf,nem,aem) +!$OMP& x,y,z,f,rpole,texfld,exf,nem,aem) !$OMP DO reduction(+:exf,nem,aem) do ii = 1, npole i = ipole(ii) @@ -339,12 +354,12 @@ subroutine exfield3 (mode,exf) yi = y(i) zi = z(i) ci = rpole(1,i) - phi = xi*exfld(1) + yi*exfld(2) + zi*exfld(3) + phi = xi*texfld(1) + yi*texfld(2) + zi*texfld(3) dix = rpole(2,i) diy = rpole(3,i) diz = rpole(4,i) - e = -f * (ci*phi + dix*exfld(1) - & + diy*exfld(2) + diz*exfld(3)) + e = -f * (ci*phi + dix*texfld(1) + & + diy*texfld(2) + diz*texfld(3)) exf = exf + e nem = nem + 1 aem(i) = aem(i) + e diff --git a/source/extfld.f b/source/extfld.f index d4d99a3bc..b02acf58b 100644 --- a/source/extfld.f +++ b/source/extfld.f @@ -12,13 +12,19 @@ c ################################################################# c c -c exfld components of applied external electric field -c use_exfld flag to include applied external electric field +c exfreq frequency of applied external field in gigahertz +c exfld components of applied external electric field +c texfld components of time dependent applied electric field +c use_exfld flag to include applied external electric field +c use_exfreq flag to oscillate applied external field c c module extfld implicit none + real*8 exfreq real*8 exfld(3) + real*8 texfld(3) logical use_exfld + logical use_exfreq save end diff --git a/source/final.f b/source/final.f index 9be300d30..17252eeef 100644 --- a/source/final.f +++ b/source/final.f @@ -88,6 +88,7 @@ subroutine final use merck use molcul use moldyn + use moment use mpole use mrecip use mutant @@ -97,6 +98,7 @@ subroutine final use opbend use opdist use orbits + use output use params use paths use pbstuf @@ -814,6 +816,10 @@ subroutine final if (allocated(aslow)) deallocate (aslow) if (allocated(afast)) deallocate (afast) c +c deallocation of global arrays from module moment +c + if (allocated(momuse)) deallocate (momuse) +c c deallocation of global arrays from module mpole c if (allocated(ipole)) deallocate (ipole) @@ -890,6 +896,12 @@ subroutine final if (allocated(worb)) deallocate (worb) if (allocated(emorb)) deallocate (emorb) c +c deallocation of global arrays from module output +c + if (allocated(ionly)) deallocate (ionly) + if (allocated(ionlyinv)) deallocate (ionlyinv) + if (allocated(print3n)) deallocate (print3n) +c c deallocation of global arrays from module params c if (allocated(prmline)) deallocate (prmline) diff --git a/source/induce.f b/source/induce.f index 7dba0e394..590ff517a 100644 --- a/source/induce.f +++ b/source/induce.f @@ -219,8 +219,8 @@ subroutine induce0a do ii = 1, npole i = ipole(ii) do j = 1, 3 - field(j,i) = field(j,i) + exfld(j) - fieldp(j,i) = fieldp(j,i) + exfld(j) + field(j,i) = field(j,i) + texfld(j) + fieldp(j,i) = fieldp(j,i) + texfld(j) end do end do end if @@ -259,7 +259,7 @@ subroutine induce0a end do end if end do - +c c get induced dipoles via the OPT extrapolation method c if (poltyp .eq. 'OPT') then @@ -4640,10 +4640,10 @@ subroutine induce0c do ii = 1, npole i = ipole(ii) do j = 1, 3 - field(j,i) = field(j,i) + exfld(j) - fieldp(j,i) = fieldp(j,i) + exfld(j) - fields(j,i) = fields(j,i) + exfld(j) - fieldps(j,i) = fieldps(j,i) + exfld(j) + field(j,i) = field(j,i) + texfld(j) + fieldp(j,i) = fieldp(j,i) + texfld(j) + fields(j,i) = fields(j,i) + texfld(j) + fieldps(j,i) = fieldps(j,i) + texfld(j) end do end do end if @@ -6011,10 +6011,10 @@ subroutine induce0d do ii = 1, npole i = ipole(ii) do j = 1, 3 - field(j,i) = field(j,i) + exfld(j) - fieldp(j,i) = fieldp(j,i) + exfld(j) - fields(j,i) = fields(j,i) + exfld(j) - fieldps(j,i) = fieldps(j,i) + exfld(j) + field(j,i) = field(j,i) + texfld(j) + fieldp(j,i) = fieldp(j,i) + texfld(j) + fields(j,i) = fields(j,i) + texfld(j) + fieldps(j,i) = fieldps(j,i) + texfld(j) end do end do end if diff --git a/source/initprm.f b/source/initprm.f index 0634201b9..c3aa9d212 100644 --- a/source/initprm.f +++ b/source/initprm.f @@ -380,9 +380,12 @@ subroutine initprm neutnbr = .false. neutcut = .false. use_exfld = .false. + use_exfreq = .false. do i = 1, 3 exfld(i) = 0.0d0 + texfld(i) = 0.0d0 end do + exfreq = 0.0d0 c c set default control parameters for atomic multipole terms c diff --git a/source/kmpole.f b/source/kmpole.f index 46309e668..c43b80155 100644 --- a/source/kmpole.f +++ b/source/kmpole.f @@ -705,6 +705,10 @@ subroutine kmpole if (use_mpole .and. .not.use_polar .and. .not.use_chgtrn) & call chkpole c +c set exclusion in computing moment of system +c + call msystem +c c turn off atomic multipole potentials if not used c if (npole .eq. 0) use_mpole = .false. diff --git a/source/mdinit.f b/source/mdinit.f index 328b6d768..efdb0f955 100644 --- a/source/mdinit.f +++ b/source/mdinit.f @@ -34,6 +34,7 @@ subroutine mdinit (dt) use moldyn use mpole use output + use polar use potent use rgddyn use rigid @@ -74,6 +75,13 @@ subroutine mdinit (dt) velsave = .false. frcsave = .false. uindsave = .false. + ustcsave = .false. + uchgsave = .false. + usyssave = .false. + vsyssave = .false. + udirsave = .false. + defsave = .false. + tefsave = .false. friction = 0.5d0 if (use_solv) friction = 91.0d0 use_sdarea = .false. @@ -125,8 +133,22 @@ subroutine mdinit (dt) velsave = .true. else if (keyword(1:11) .eq. 'SAVE-FORCE ') then frcsave = .true. - else if (keyword(1:13) .eq. 'SAVE-INDUCED ') then + else if (keyword(1:13) .eq. 'SAVE-UINDUCE ') then uindsave = .true. + else if (keyword(1:13) .eq. 'SAVE-USTATIC ') then + ustcsave = .true. + else if (keyword(1:13) .eq. 'SAVE-UCHARGE ') then + uchgsave = .true. + else if (keyword(1:13) .eq. 'SAVE-USYSTEM ') then + usyssave = .true. + else if (keyword(1:13) .eq. 'SAVE-VSYSTEM ') then + vsyssave = .true. + else if (keyword(1:13) .eq. 'SAVE-UDIRECT ') then + udirsave = .true. + else if (keyword(1:13) .eq. 'SAVE-DEFIELD ') then + defsave = .true. + else if (keyword(1:13) .eq. 'SAVE-TEFIELD ') then + tefsave = .true. else if (keyword(1:9) .eq. 'FRICTION ') then read (string,*,err=10,end=10) friction else if (keyword(1:17) .eq. 'FRICTION-SCALING ') then @@ -174,6 +196,14 @@ subroutine mdinit (dt) 10 continue end do c +c check for use of save-only keyword +c + call saveonly +c +c get unique atom types +c + call uniquetyp +c c check for use of induced dipole prediction methods c if (use_polar) call predict diff --git a/source/mdsave.f b/source/mdsave.f index f830e309f..eb3895898 100644 --- a/source/mdsave.f +++ b/source/mdsave.f @@ -23,11 +23,13 @@ subroutine mdsave (istep,dt,epot,eksum) use bound use boxes use couple + use extfld use files use group use inform use iounit use mdstuf + use moment use mpole use output use polar @@ -35,17 +37,24 @@ subroutine mdsave (istep,dt,epot,eksum) use rgddyn use socket use titles + use uatom + use units implicit none integer i,j,lext integer istep integer ixyz,iind - integer ivel,ifrc + integer ivel,ifrc,istc integer iend,isave integer freeunit integer trimtext integer modsave real*8 dt,pico real*8 epot,eksum + real*8 xustc,yustc,zustc + real*8 xuind,yuind,zuind + real*8 xuchg,yuchg,zuchg + real*8 xf,yf,zf + real*8 xm,ym,zm logical exist,first character*7 ext character*240 endfile @@ -53,6 +62,7 @@ subroutine mdsave (istep,dt,epot,eksum) character*240 velfile character*240 frcfile character*240 indfile + character*240 stcfile c c c send data via external socket communication if desired @@ -123,52 +133,147 @@ subroutine mdsave (istep,dt,epot,eksum) end if end if c +c print the external electric field values +c + if (use_exfld) then + xf = texfld(1) * elefield + yf = texfld(2) * elefield + zf = texfld(3) * elefield + if (digits .le. 6) then + write (iout,170) xf,yf,zf + 170 format (' External Field',7x,3f14.6) + else if (digits .le. 8) then + write (iout,180) xf,yf,zf + 180 format (' External Field',7x,3f16.8) + else + write (iout,190) xf,yf,zf + 190 format (' External Field',7x,3f18.10) + end if + end if +c c move stray molecules into periodic box if desired c if (use_wrap) call bounds c -c save coordinates to archive or numbered structure file +c compute center of mass if saving dipole moment c - ixyz = freeunit () - if (cyclesave) then - xyzfile = filename(1:leng)//'.'//ext(1:lext) - call version (xyzfile,'new') - open (unit=ixyz,file=xyzfile,status='new') - call prtxyz (ixyz) - else if (dcdsave) then - xyzfile = filename(1:leng) - call suffix (xyzfile,'dcd','old') - inquire (file=xyzfile,exist=exist) - if (exist) then - first = .false. - open (unit=ixyz,file=xyzfile,form='unformatted', - & status='old',position='append') + if (usyssave .or. uchgsave) call compcent(xm,ym,zm) +c +c compute dipole moment of system +c + if (usyssave) then + call dmoments (xm,ym,zm,xustc,yustc,zustc,xuind,yuind,zuind, + & xuchg,yuchg,zuchg) + if (digits .le. 6) then + write (iout,200) xuchg,yuchg,zuchg + 200 format (' System Charge Dipole',1x,3f20.6) + write (iout,210) xustc,yustc,zustc + 210 format (' System Static Dipole',1x,3f20.6) + else if (digits .le. 8) then + write (iout,220) xuchg,yuchg,zuchg + 220 format (' System Charge Dipole',1x,3f22.8) + write (iout,230) xustc,yustc,zustc + 230 format (' System Static Dipole',1x,3f22.8) else - first = .true. - open (unit=ixyz,file=xyzfile,form='unformatted', - & status='new') + write (iout,240) xuchg,yuchg,zuchg + 240 format (' System Charge Dipole',1x,3f24.10) + write (iout,250) xustc,yustc,zustc + 250 format (' System Static Dipole',1x,3f24.10) end if - call prtdcd (ixyz,first) - else - xyzfile = filename(1:leng) - call suffix (xyzfile,'arc','old') - inquire (file=xyzfile,exist=exist) - if (exist) then - call openend (ixyz,xyzfile) - else + if (use_polar) then + if (digits .le. 6) then + write (iout,260) xuind,yuind,zuind + 260 format (' System Induced Dipole',3f20.6) + else if (digits .le. 8) then + write (iout,270) xuind,yuind,zuind + 270 format (' System Induced Dipole',3f22.8) + else + write (iout,280) xuind,yuind,zuind + 280 format (' System Induced Dipole',3f24.10) + end if + end if + write (iout,290) + 290 format (' Charge Dipole by Atom Type:', + & /,' Type',11x,'X-UCharge',11x,'Y-UCharge',11x,'Z-UCharge') + do i = 1, nunique + write (iout,300) utype(i),utv3(1,i),utv3(2,i),utv3(3,i) + 300 format (i5,3f20.6) + end do + write (iout,310) + 310 format (' Static Dipole by Atom Type:', + & /,' Type',11x,'X-UStatic',11x,'Y-UStatic',11x,'Z-UStatic') + do i = 1, nunique + write (iout,320) utype(i),utv1(1,i),utv1(2,i),utv1(3,i) + 320 format (i5,3f20.6) + end do + if (use_polar) then + write (iout,330) + 330 format (' Induced Dipole by Atom Type:', + & /,' Type',11x,'X-UInduce',11x,'Y-UInduce',11x,'Z-UInduce') + do i = 1, nunique + write (iout,340) utype(i),utv2(1,i),utv2(2,i),utv2(3,i) + 340 format (i5,3f20.6) + end do + end if + end if +c +c compute velocity of unique atom types in the system +c + if (vsyssave) then + call velunique + write (iout,350) + 350 format (' Velocity by Atom Type:', + & /,' Type',10x,'X-Velocity',10x,'Y-Velocity',10x,'Z-Velocity') + do i = 1, nunique + write (iout,360) utype(i),utv1(1,i),utv1(2,i),utv1(3,i) + 360 format (i5,3f20.6) + end do + end if +c +c save coordinates to archive or numbered structure file +c + write (iout,370) isave + 370 format (' Frame Number',13x,i10) + if (coordsave) then + ixyz = freeunit () + if (cyclesave) then + xyzfile = filename(1:leng)//'.'//ext(1:lext) + call version (xyzfile,'new') open (unit=ixyz,file=xyzfile,status='new') + call prtxyz (ixyz) + else if (dcdsave) then + xyzfile = filename(1:leng) + call suffix (xyzfile,'dcd','old') + inquire (file=xyzfile,exist=exist) + if (exist) then + first = .false. + open (unit=ixyz,file=xyzfile,form='unformatted', + & status='old',position='append') + else + first = .true. + open (unit=ixyz,file=xyzfile,form='unformatted', + & status='new') + end if + call prtdcd (ixyz,first) + else + xyzfile = filename(1:leng) + call suffix (xyzfile,'arc','old') + inquire (file=xyzfile,exist=exist) + if (exist) then + call openend (ixyz,xyzfile) + else + open (unit=ixyz,file=xyzfile,status='new') + end if + call prtxyz (ixyz) end if - call prtxyz (ixyz) + close (unit=ixyz) + write (iout,380) xyzfile(1:trimtext(xyzfile)) + 380 format (' Coordinate File',13x,a) end if - close (unit=ixyz) - write (iout,170) isave - 170 format (' Frame Number',13x,i10) - write (iout,180) xyzfile(1:trimtext(xyzfile)) - 180 format (' Coordinate File',13x,a) c c update the information needed to restart the trajectory c - call prtdyn + if (dynsave) call prtdyn c c save the velocity vector components at the current step c @@ -202,22 +307,22 @@ subroutine mdsave (istep,dt,epot,eksum) end if end if if (integrate .eq. 'RIGIDBODY') then - write (ivel,190) ngrp,title(1:ltitle) - 190 format (i6,2x,a) + write (ivel,390) ngrp,title(1:ltitle) + 390 format (i6,2x,a) do i = 1, ngrp - write (ivel,200) i,(vcm(j,i),j=1,3) - 200 format (i6,3x,d13.6,3x,d13.6,3x,d13.6) - write (ivel,210) i,(wcm(j,i),j=1,3) - 210 format (i6,3x,d13.6,3x,d13.6,3x,d13.6) + write (ivel,400) i,(vcm(j,i),j=1,3) + 400 format (i6,3x,d13.6,3x,d13.6,3x,d13.6) + write (ivel,410) i,(wcm(j,i),j=1,3) + 410 format (i6,3x,d13.6,3x,d13.6,3x,d13.6) end do else if (dcdsave) then - call prtdcdv (ivel,first) + call prtdcdv3 (ivel,first,'VEL') else - call prtvel (ivel) + call prtvec3 (ivel,'VEL') end if close (unit=ivel) - write (iout,240) velfile(1:trimtext(velfile)) - 240 format (' Velocity File',15x,a) + write (iout,420) velfile(1:trimtext(velfile)) + 420 format (' Velocity File',15x,a) end if c c save the force vector components for the current step; not @@ -245,7 +350,7 @@ subroutine mdsave (istep,dt,epot,eksum) open (unit=ifrc,file=frcfile,form='unformatted', & status='new') end if - call prtdcdf (ifrc,first) + call prtdcdv3 (ifrc,first,'FRC') else frcfile = filename(1:leng) call suffix (frcfile,'frc','old') @@ -255,11 +360,11 @@ subroutine mdsave (istep,dt,epot,eksum) else open (unit=ifrc,file=frcfile,status='new') end if - call prtfrc (ifrc) + call prtvec3 (ifrc,'FRC') end if close (unit=ifrc) - write (iout,270) frcfile(1:trimtext(frcfile)) - 270 format (' Force Vector File',11x,a) + write (iout,430) frcfile(1:trimtext(frcfile)) + 430 format (' Force Vector File',11x,a) end if c c save the induced dipole components for the current step @@ -267,12 +372,13 @@ subroutine mdsave (istep,dt,epot,eksum) if (uindsave .and. use_polar) then iind = freeunit () if (cyclesave) then - indfile = filename(1:leng)//'.'//ext(1:lext)//'u' + indfile = filename(1:leng)//'.'//ext(1:lext)//'ui' call version (indfile,'new') open (unit=iind,file=indfile,status='new') + call prtvec3 (iind,'UIN') else if (dcdsave) then indfile = filename(1:leng) - call suffix (indfile,'dcdu','old') + call suffix (indfile,'dcdui','old') inquire (file=indfile,exist=exist) if (exist) then first = .false. @@ -283,7 +389,7 @@ subroutine mdsave (istep,dt,epot,eksum) open (unit=iind,file=indfile,form='unformatted', & status='new') end if - call prtdcdu (iind,first) + call prtdcdv3 (iind,first,'UIN') else indfile = filename(1:leng) call suffix (indfile,'uind','old') @@ -293,11 +399,206 @@ subroutine mdsave (istep,dt,epot,eksum) else open (unit=iind,file=indfile,status='new') end if - call prtuind (iind) + call prtvec3 (iind,'UIN') end if close (unit=iind) - write (iout,300) indfile(1:trimtext(indfile)) - 300 format (' Induced Dipole File',9x,a) + write (iout,440) indfile(1:trimtext(indfile)) + 440 format (' Induced Dipole File',9x,a) + end if +c +c save the static dipole components for the current step +c + if (ustcsave) then + istc = freeunit () + if (cyclesave) then + stcfile = filename(1:leng)//'.'//ext(1:lext)//'us' + call version (stcfile,'new') + open (unit=istc,file=stcfile,status='new') + call prtvec3 (istc,'UST') + else if (dcdsave) then + stcfile = filename(1:leng) + call suffix (stcfile,'dcdus','old') + inquire (file=stcfile,exist=exist) + if (exist) then + first = .false. + open (unit=istc,file=stcfile,form='unformatted', + & status='old',position='append') + else + first = .true. + open (unit=istc,file=stcfile,form='unformatted', + & status='new') + end if + call prtdcdv3 (istc,first,'UST') + else + stcfile = filename(1:leng) + call suffix (stcfile,'ustc','old') + inquire (file=stcfile,exist=exist) + if (exist) then + call openend (istc,stcfile) + else + open (unit=istc,file=stcfile,status='new') + end if + call prtvec3 (istc,'UST') + end if + close (unit=istc) + write (iout,450) stcfile(1:trimtext(stcfile)) + 450 format (' Static Dipole File',10x,a) + end if +c +c save the charge dipole components for the current step +c + if (uchgsave) then + istc = freeunit () + if (cyclesave) then + stcfile = filename(1:leng)//'.'//ext(1:lext)//'uc' + call version (stcfile,'new') + open (unit=istc,file=stcfile,status='new') + call prtvec3 (istc,'UCH') + else if (dcdsave) then + stcfile = filename(1:leng) + call suffix (stcfile,'dcduc','old') + inquire (file=stcfile,exist=exist) + if (exist) then + first = .false. + open (unit=istc,file=stcfile,form='unformatted', + & status='old',position='append') + else + first = .true. + open (unit=istc,file=stcfile,form='unformatted', + & status='new') + end if + call prtdcdv3 (istc,first,'UCH') + else + stcfile = filename(1:leng) + call suffix (stcfile,'uchg','old') + inquire (file=stcfile,exist=exist) + if (exist) then + call openend (istc,stcfile) + else + open (unit=istc,file=stcfile,status='new') + end if + call prtvec3 (istc,'UCH') + end if + close (unit=istc) + write (iout,460) stcfile(1:trimtext(stcfile)) + 460 format (' Charge Dipole File',10x,a) + end if +c +c save the direct induced dipole components for the current step +c + if (udirsave .and. use_polar) then + iind = freeunit () + if (cyclesave) then + indfile = filename(1:leng)//'.'//ext(1:lext)//'ud' + call version (indfile,'new') + open (unit=iind,file=indfile,status='new') + call prtvec3 (iind,'UDR') + else if (dcdsave) then + indfile = filename(1:leng) + call suffix (indfile,'dcdud','old') + inquire (file=indfile,exist=exist) + if (exist) then + first = .false. + open (unit=iind,file=indfile,form='unformatted', + & status='old',position='append') + else + first = .true. + open (unit=iind,file=indfile,form='unformatted', + & status='new') + end if + call prtdcdv3 (iind,first,'UDR') + else + indfile = filename(1:leng) + call suffix (indfile,'udir','old') + inquire (file=indfile,exist=exist) + if (exist) then + call openend (iind,indfile) + else + open (unit=iind,file=indfile,status='new') + end if + call prtvec3 (iind,'UDR') + end if + close (unit=iind) + write (iout,470) indfile(1:trimtext(indfile)) + 470 format (' Direct Induced Dipole File',2x,a) + end if +c +c save the direct atomic electric field for the current step +c + if (defsave .and. use_polar) then + iind = freeunit () + if (cyclesave) then + indfile = filename(1:leng)//'.'//ext(1:lext)//'de' + call version (indfile,'new') + open (unit=iind,file=indfile,status='new') + call prtvec3 (iind,'DEF') + else if (dcdsave) then + indfile = filename(1:leng) + call suffix (indfile,'dcdde','old') + inquire (file=indfile,exist=exist) + if (exist) then + first = .false. + open (unit=iind,file=indfile,form='unformatted', + & status='old',position='append') + else + first = .true. + open (unit=iind,file=indfile,form='unformatted', + & status='new') + end if + call prtdcdv3 (iind,first,'DEF') + else + indfile = filename(1:leng) + call suffix (indfile,'def','old') + inquire (file=indfile,exist=exist) + if (exist) then + call openend (iind,indfile) + else + open (unit=iind,file=indfile,status='new') + end if + call prtvec3 (iind,'DEF') + end if + close (unit=iind) + write (iout,480) indfile(1:trimtext(indfile)) + 480 format (' Direct Electric Field File',2x,a) + end if +c +c save the total atomic electric field for the current step +c + if (tefsave .and. use_polar) then + iind = freeunit () + if (cyclesave) then + indfile = filename(1:leng)//'.'//ext(1:lext)//'te' + call version (indfile,'new') + open (unit=iind,file=indfile,status='new') + call prtvec3 (iind,'TEF') + else if (dcdsave) then + indfile = filename(1:leng) + call suffix (indfile,'dcdte','old') + inquire (file=indfile,exist=exist) + if (exist) then + first = .false. + open (unit=iind,file=indfile,form='unformatted', + & status='old',position='append') + else + first = .true. + open (unit=iind,file=indfile,form='unformatted', + & status='new') + end if + call prtdcdv3 (iind,first,'TEF') + else + indfile = filename(1:leng) + call suffix (indfile,'tef','old') + inquire (file=indfile,exist=exist) + if (exist) then + call openend (iind,indfile) + else + open (unit=iind,file=indfile,status='new') + end if + call prtvec3 (iind,'TEF') + end if + close (unit=iind) + write (iout,490) indfile(1:trimtext(indfile)) + 490 format (' Total Electric Field File',3x,a) end if c c test for requested termination of the dynamics calculation @@ -314,8 +615,8 @@ subroutine mdsave (istep,dt,epot,eksum) end if end if if (exist) then - write (iout,310) - 310 format (/,' MDSAVE -- Dynamics Calculation Ending', + write (iout,500) + 500 format (/,' MDSAVE -- Dynamics Calculation Ending', & ' due to User Request') call fatal end if @@ -324,8 +625,8 @@ subroutine mdsave (istep,dt,epot,eksum) c modsave = mod(istep,iprint) if (verbose .and. modsave.ne.0) then - write (iout,320) - 320 format () + write (iout,510) + 510 format () end if return end diff --git a/source/moment.f b/source/moment.f index d6046ff26..840e80429 100644 --- a/source/moment.f +++ b/source/moment.f @@ -27,6 +27,7 @@ c zxqpl total quadrupole tensor zx-component in global frame c zyqpl total quadrupole tensor zy-component in global frame c zzqpl total quadrupole tensor zz-component in global frame +c momuse true if an atom is used for system moment calculation c c module moment @@ -37,5 +38,6 @@ module moment real*8 xxqpl,xyqpl,xzqpl real*8 yxqpl,yyqpl,yzqpl real*8 zxqpl,zyqpl,zzqpl + logical, allocatable :: momuse(:) save end diff --git a/source/moments.f b/source/moments.f index d9dfced0d..e83c163c7 100644 --- a/source/moments.f +++ b/source/moments.f @@ -131,7 +131,7 @@ subroutine moments (mode) c do i = 1, nion k = iion(i) - if (use(k)) then + if (use(k) .and. momuse(k)) then netchg = netchg + pchg(k) xdpl = xdpl + xcm(k)*pchg(k) ydpl = ydpl + ycm(k)*pchg(k) @@ -153,7 +153,7 @@ subroutine moments (mode) do i = 1, ndipole j = idpl(1,i) k = idpl(2,i) - if (use(j) .or. use(k)) then + if ((use(j).and.momuse(j)) .or. (use(k).and.momuse(k))) then xi = x(j) - x(k) yi = y(j) - y(k) zi = z(j) - z(k) @@ -205,7 +205,7 @@ subroutine moments (mode) c do i = 1, npole k = ipole(i) - if (use(k)) then + if (use(k) .and. momuse(k)) then netchg = netchg + rpole(1,k) xdpl = xdpl + xcm(k)*rpole(1,k) + rpole(2,k) ydpl = ydpl + ycm(k)*rpole(1,k) + rpole(3,k) @@ -254,7 +254,7 @@ subroutine moments (mode) c do i = 1, npole k = ipole(i) - if (use(k)) then + if (use(k) .and. momuse(k)) then xxqpl = xxqpl + 3.0d0*rpole(5,k) xyqpl = xyqpl + 3.0d0*rpole(6,k) xzqpl = xzqpl + 3.0d0*rpole(7,k) @@ -309,3 +309,220 @@ subroutine moments (mode) call jacobi (3,a,netqpl,b) return end +c +c +c ############################################################## +c ## ## +c ## subroutine dmoments -- total electric dipole moments ## +c ## ## +c ############################################################## +c +c +c "dmoments" computes the total dipole moments over all atoms and +c the unique atom type dipole moments; +c called in mdsave, it is assumed bound is called +c +c + subroutine dmoments (xm,ym,zm,xustc,yustc,zustc,xuind,yuind,zuind, + & xuchg,yuchg,zuchg) + use atoms + use dipole + use charge + use moment + use mpole + use polar + use potent + use solpot + use uatom + use units + implicit none + integer i,j,k + integer ii + integer ut + real*8 xi,yi,zi,ri + real*8 xm,ym,zm + real*8 xu,yu,zu + real*8 xc,yc,zc + real*8 xustc,yustc,zustc + real*8 xuind,yuind,zuind + real*8 xuchg,yuchg,zuchg + real*8 xbnd,ybnd,zbnd +c +c +c zero out dipole moments +c + xustc = 0.0d0 + yustc = 0.0d0 + zustc = 0.0d0 + xuind = 0.0d0 + yuind = 0.0d0 + zuind = 0.0d0 + xuchg = 0.0d0 + yuchg = 0.0d0 + zuchg = 0.0d0 + do i = 1, nunique + do j = 1, 3 + utv1(j,i) = 0.0d0 + utv2(j,i) = 0.0d0 + utv3(j,i) = 0.0d0 + end do + end do +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(nion,iion,x,y,z,pchg,momuse,xuchg,yuchg,zuchg,utv3,type, +!$OMP& utypeinv) +!$OMP DO reduction(+:xuchg,yuchg,zuchg,utv3) +c +c compute the static dipole moment due to partial charges +c + do ii = 1, nion + i = iion(ii) + if (momuse(i)) then + xc = (x(i)-xm)*pchg(i) + yc = (y(i)-ym)*pchg(i) + zc = (z(i)-zm)*pchg(i) + xuchg = xuchg + xc + yuchg = yuchg + yc + zuchg = zuchg + zc + ut = utypeinv(type(i)) + utv3(1,ut) = utv3(1,ut) + xc + utv3(2,ut) = utv3(2,ut) + yc + utv3(3,ut) = utv3(3,ut) + zc + end if + end do +!$OMP END DO +!$OMP END PARALLEL +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(ndipole,x,y,z,idpl,bdpl,momuse,xustc,yustc,zustc, +!$OMP& utv1,type,utypeinv) +!$OMP DO reduction(+:xustc,yustc,zustc,utv1) +c +c compute the static dipole moment due to bond dipoles +c + do i = 1, ndipole + j = idpl(1,i) + k = idpl(2,i) + if (momuse(j) .or. momuse(k)) then + xi = x(j) - x(k) + yi = y(j) - y(k) + zi = z(j) - z(k) + ri = sqrt(xi*xi + yi*yi + zi*zi) + xbnd = bdpl(i) * (xi/ri) / debye + ybnd = bdpl(i) * (yi/ri) / debye + zbnd = bdpl(i) * (zi/ri) / debye + xustc = xustc + xbnd + yustc = yustc + ybnd + zustc = zustc + zbnd + ut = utypeinv(type(j)) + utv1(1,ut) = utv1(1,ut) + xbnd + utv1(2,ut) = utv1(2,ut) + ybnd + utv1(3,ut) = utv1(3,ut) + zbnd + end if + end do +!$OMP END DO +!$OMP END PARALLEL +c +c OpenMP directives for the major loop structure +c +!$OMP PARALLEL default(private) +!$OMP& shared(npole,x,y,z,ipole,rpole,uind,uinds,use_polar,solvtyp, +!$OMP& momuse,xm,ym,zm,xustc,yustc,zustc,xuind,yuind,zuind, +!$OMP& xuchg,yuchg,zuchg,utv1,utv2,utv3,type,utypeinv) +!$OMP DO reduction(+:xustc,yustc,zustc,xuchg,yuchg,zuchg,utv1,utv3) +c +c compute the static dipole moment due to atomic multipoles +c + do ii = 1, npole + i = ipole(ii) + if (momuse(i)) then + xu = rpole(2,i) + yu = rpole(3,i) + zu = rpole(4,i) + xc = (x(i)-xm)*rpole(1,i) + yc = (y(i)-ym)*rpole(1,i) + zc = (z(i)-zm)*rpole(1,i) + xustc = xustc + xu + yustc = yustc + yu + zustc = zustc + zu + xuchg = xuchg + xc + yuchg = yuchg + yc + zuchg = zuchg + zc + ut = utypeinv(type(i)) + utv1(1,ut) = utv1(1,ut) + xu + utv1(2,ut) = utv1(2,ut) + yu + utv1(3,ut) = utv1(3,ut) + zu + utv3(1,ut) = utv3(1,ut) + xc + utv3(2,ut) = utv3(2,ut) + yc + utv3(3,ut) = utv3(3,ut) + zc + end if + end do +!$OMP END DO +c +c compute the induced dipole moment +c + if (use_polar) then + if (solvtyp.eq.'GK' .or. solvtyp.eq.'PB') then +!$OMP DO reduction(+:xuind,yuind,zuind,utv2) + do ii = 1, npole + i = ipole(ii) + if (momuse(i)) then + xu = uinds(1,i) + yu = uinds(2,i) + zu = uinds(3,i) + xuind = xuind + xu + yuind = yuind + yu + zuind = zuind + zu + ut = utypeinv(type(i)) + utv2(1,ut) = utv2(1,ut) + xu + utv2(2,ut) = utv2(2,ut) + yu + utv2(3,ut) = utv2(3,ut) + zu + end if + end do +!$OMP END DO + else +!$OMP DO reduction(+:xuind,yuind,zuind,utv2) + do ii = 1, npole + i = ipole(ii) + if (momuse(i)) then + xu = uind(1,i) + yu = uind(2,i) + zu = uind(3,i) + xuind = xuind + xu + yuind = yuind + yu + zuind = zuind + zu + ut = utypeinv(type(i)) + utv2(1,ut) = utv2(1,ut) + xu + utv2(2,ut) = utv2(2,ut) + yu + utv2(3,ut) = utv2(3,ut) + zu + end if + end do +!$OMP END DO + end if + end if +!$OMP END PARALLEL +c +c convert dipole to Debye +c + xustc = xustc * debye + yustc = yustc * debye + zustc = zustc * debye + xuind = xuind * debye + yuind = yuind * debye + zuind = zuind * debye + xuchg = xuchg * debye + yuchg = yuchg * debye + zuchg = zuchg * debye + do i = 1, nunique + do j = 1, 3 + utv1(j,i) = utv1(j,i) * debye + utv2(j,i) = utv2(j,i) * debye + utv3(j,i) = utv3(j,i) * debye + end do + end do + return + end diff --git a/source/optinit.f b/source/optinit.f index b930db14e..3f52c88ef 100644 --- a/source/optinit.f +++ b/source/optinit.f @@ -53,7 +53,7 @@ subroutine optinit use_wrap = .true. else if (keyword(1:11) .eq. 'SAVE-FORCE ') then frcsave = .true. - else if (keyword(1:13) .eq. 'SAVE-INDUCED ') then + else if (keyword(1:13) .eq. 'SAVE-UINDUCE ') then uindsave = .true. end if 10 continue diff --git a/source/optsave.f b/source/optsave.f index d97129d81..5679b139d 100644 --- a/source/optsave.f +++ b/source/optsave.f @@ -175,7 +175,7 @@ subroutine optsave (ncycle,f,xx) open (unit=ifrc,file=frcfile,form='unformatted', & status='new') end if - call prtdcdf (ifrc,first) + call prtdcdv3 (ifrc,first,'FRC') else if (arcsave) then frcfile = filename(1:leng) call suffix (frcfile,'frc','old') @@ -185,12 +185,12 @@ subroutine optsave (ncycle,f,xx) else open (unit=ifrc,file=frcfile,status='new') end if - call prtfrc (ifrc) + call prtvec3 (ifrc,'FRC') else frcfile = filename(1:leng)//'.'//ext(1:lext)//'f' call version (frcfile,'new') open (unit=ifrc,file=frcfile,status='new') - call prtfrc (ifrc) + call prtvec3 (ifrc,'FRC') end if else frcfile = filename(1:leng) @@ -203,7 +203,7 @@ subroutine optsave (ncycle,f,xx) open (unit=ifrc,file=frcfile,status='new') end if rewind (unit=ifrc) - call prtfrc (ifrc) + call prtvec3 (ifrc,'FRC') end if close (unit=ifrc) write (iout,10) frcfile(1:trimtext(frcfile)) @@ -228,7 +228,7 @@ subroutine optsave (ncycle,f,xx) open (unit=iind,file=indfile,form='unformatted', & status='new') end if - call prtdcdu (iind,first) + call prtdcdv3 (iind,first,'UIN') else if (arcsave) then indfile = filename(1:leng) call suffix (indfile,'uind','old') @@ -238,12 +238,12 @@ subroutine optsave (ncycle,f,xx) else open (unit=iind,file=indfile,status='new') end if - call prtuind (iind) + call prtvec3 (iind,'UIN') else indfile = filename(1:leng)//'.'//ext(1:lext)//'u' call version (indfile,'new') open (unit=iind,file=indfile,status='new') - call prtuind (iind) + call prtvec3 (iind,'UIN') end if else indfile = filename(1:leng) @@ -256,7 +256,7 @@ subroutine optsave (ncycle,f,xx) open (unit=iind,file=indfile,status='new') end if rewind (unit=iind) - call prtuind (iind) + call prtvec3 (iind,'UIN') end if close (unit=iind) write (iout,20) indfile(1:trimtext(indfile)) diff --git a/source/output.f b/source/output.f index 9d36de979..39151c7e4 100644 --- a/source/output.f +++ b/source/output.f @@ -12,31 +12,59 @@ c ################################################################ c c +c nonly total number of save sites in the system +c ionly number of the atom for each save site +c ionlyinv inverse map of ionly +c print3n array for printing atomic 3D vectors c archive logical flag for coordinates in Tinker XYZ format c binary logical flag for coordinates in DCD binary format c noversion logical flag governing use of filename versions c overwrite logical flag to overwrite intermediate files inplace +c coordsave logical flag to save coordinates +c dynsave logical flag to save dynamics (.dyn) file +c onlysave logical flag to only save certain coordinates c arcsave logical flag to save coordinates in Tinker XYZ format c dcdsave logical flag to save coordinates in DCD binary format c cyclesave logical flag to mark use of numbered cycle files c velsave logical flag to save velocity vector components c frcsave logical flag to save force vector components c uindsave logical flag to save induced atomic dipoles +c ustcsave logical flag to save static atomic dipoles +c uchgsave logical flag to save charge atomic dipoles +c usyssave logical flag to save unique atom type dipole moment +c vsyssave logical flag to save unique atom type velocity +c udirsave logical flag to save the direct induced atomic dipoles +c defsave logical flag to save the direct electric field +c tefsave logical flag to save the total electric field c coordtype selects Cartesian, internal, rigid body or none c c module output implicit none + integer nonly + integer, allocatable :: ionly(:) + integer, allocatable :: ionlyinv(:) + real*8, allocatable :: print3n(:,:) logical archive logical binary logical noversion logical overwrite + logical coordsave + logical dynsave logical cyclesave + logical onlysave logical arcsave logical dcdsave logical velsave logical frcsave logical uindsave + logical ustcsave + logical uchgsave + logical usyssave + logical vsyssave + logical udirsave + logical defsave + logical tefsave character*9 coordtype save end diff --git a/source/prmkey.f b/source/prmkey.f index dd81b90b0..2d77b91db 100644 --- a/source/prmkey.f +++ b/source/prmkey.f @@ -25,6 +25,7 @@ subroutine prmkey (text) use expol use extfld use fields + use math use mplpot use polpot use potent @@ -403,7 +404,12 @@ subroutine prmkey (text) use_exfld = .true. do i = 1, 3 exfld(i) = exfld(i) / elefield + texfld(i) = exfld(i) end do + else if (keyword(1:15) .eq. 'EXFLD-FREQ ') then + read (string,*,err=10,end=10) exfreq + use_exfreq = .true. + exfreq = 2.0d0 * pi * 0.001d0 * exfreq c c set control parameters for atomic multipole potentials c diff --git a/source/prtarc.f b/source/prtarc.f index 970d56057..b98d7e7fa 100644 --- a/source/prtarc.f +++ b/source/prtarc.f @@ -240,10 +240,12 @@ subroutine prtarcb (idcd,first) write (idcd) nuse end if c -c append the lattice values based on header flag value +c append the lattice values based on header flag value; +c using angle values is NAMD style, cosine values is CHARMM c if (use_bounds) then - write (idcd) xbox,gamma_cos,ybox,beta_cos,alpha_cos,zbox +c write (idcd) xbox,gamma_cos,ybox,beta_cos,alpha_cos,zbox + write (idcd) xbox,gamma,ybox,beta,alpha,zbox end if c c remove unused atoms from the coordinates to be output diff --git a/source/prtfrc.f b/source/prtfrc.f deleted file mode 100644 index 18c08fe9e..000000000 --- a/source/prtfrc.f +++ /dev/null @@ -1,230 +0,0 @@ -c -c -c ################################################### -c ## COPYRIGHT (C) 2023 by Jay William Ponder ## -c ## All Rights Reserved ## -c ################################################### -c -c ############################################################## -c ## ## -c ## subroutine prtfrc -- output of atom force components ## -c ## ## -c ############################################################## -c -c -c "prtfrc" writes out a set of atom force components to an -c external file in Tinker XYZ format -c -c - subroutine prtfrc (ifrc) - use atomid - use atoms - use bound - use boxes - use couple - use deriv - use files - use inform - use titles - implicit none - integer i,j,k,ifrc - integer size,crdsiz - real*8 crdmin,crdmax - logical opened - character*2 atmc - character*2 crdc - character*2 digc - character*25 fstr - character*240 frcfile -c -c -c open the output unit if not already done -c - inquire (unit=ifrc,opened=opened) - if (.not. opened) then - frcfile = filename(1:leng)//'.frc' - call version (frcfile,'new') - open (unit=ifrc,file=frcfile,status='new') - end if -c -c check for large systems needing extended formatting -c - atmc = 'i6' - if (n .ge. 100000) atmc = 'i7' - if (n .ge. 1000000) atmc = 'i8' - crdmin = 0.0d0 - crdmax = 0.0d0 - do i = 1, n - crdmin = min(crdmin,x(i),y(i),z(i)) - crdmax = max(crdmax,x(i),y(i),z(i)) - end do - crdsiz = 6 - if (crdmin .le. -1000.0d0) crdsiz = 7 - if (crdmax .ge. 10000.0d0) crdsiz = 7 - if (crdmin .le. -10000.0d0) crdsiz = 8 - if (crdmax .ge. 100000.0d0) crdsiz = 8 - crdsiz = crdsiz + max(6,digits) - size = 0 - call numeral (crdsiz,crdc,size) - if (digits .le. 6) then - digc = '6 ' - else if (digits .le. 8) then - digc = '8' - else - digc = '10' - end if -c -c write out the number of atoms and the title -c - if (ltitle .eq. 0) then - fstr = '('//atmc//')' - write (ifrc,fstr(1:4)) n - else - fstr = '('//atmc//',2x,a)' - write (ifrc,fstr(1:9)) n,title(1:ltitle) - end if -c -c write out the periodic cell lengths and angles -c - if (use_bounds) then - fstr = '(1x,6f'//crdc//'.'//digc//')' - write (ifrc,fstr) xbox,ybox,zbox,alpha,beta,gamma - end if -c -c write out the atom force components for each atom -c - fstr = '('//atmc//',2x,a3,3f'//crdc// - & '.'//digc//',i6,8'//atmc//')' - do i = 1, n - k = n12(i) - if (k .eq. 0) then - write (ifrc,fstr) i,name(i),(-desum(j,i),j=1,3),type(i) - else - write (ifrc,fstr) i,name(i),(-desum(j,i),j=1,3),type(i), - & (i12(j,i),j=1,k) - end if - end do -c -c close the output unit if opened by this routine -c - if (.not. opened) close (unit=ifrc) - return - end -c -c -c ############################################################## -c ## ## -c ## subroutine prtdcdf -- output of DCD force components ## -c ## ## -c ############################################################## -c -c -c "prtdcdf" writes out a set of atomic force components to -c a file in CHARMM DCD binary format compatible with the VMD -c visualization software and other packages -c -c note the format used is based on the "dcdplugin.c" code from -c the NAMD and VMD programs, and tutorial 4.1 from the software -c package GENESIS: Generalized-Ensemble Simulation System -c -c variables and parameters: -c -c header type of data (CORD=coordinates, VELD=velocities) -c nframe number of frames stored in the DCD file -c nprev number of previous integration steps -c ncrdsav frequency in steps for saving coordinate frames -c nstep number of integration steps in the total run -c nvelsav frequency of coordinate saves with velocity data -c ndfree number of degrees of freedom for the system -c nfixat number of fixed atoms for the system -c usebox flag for periodic boundaries (1=true, 0=false) -c use4d flag for 4D trajectory (1=true, 0=false) -c usefq flag for fluctuating charges (1=true, 0=false) -c merged result of merge without checks (1=true, 0=false) -c vcharmm version of CHARMM software for compatibility -c -c in general a value of zero for any of the above indicates that -c the particular feature is unused -c -c - subroutine prtdcdf (idcd,first) - use atoms - use bound - use boxes - use deriv - use files - use titles - implicit none - integer i,idcd - integer zero,one - integer nframe,nprev - integer ncrdsav,nstep - integer nvelsav,ndfree - integer nfixat,usebox - integer use4d,usefq - integer merged,vcharmm - integer ntitle - real*4 tdelta - logical opened,first - character*4 header - character*240 dcdfile -c -c -c open the output unit if not already done -c - inquire (unit=idcd,opened=opened) - if (.not. opened) then - dcdfile = filename(1:leng)//'.dcdf' - call version (dcdfile,'new') - open (unit=idcd,file=dcdfile,form='unformatted',status='new') - end if -c -c write header info along with title and number of atoms -c - if (first) then - first = .false. - zero = 0 - one = 1 - header = 'CORD' - nframe = zero - nprev = zero - ncrdsav = one - nstep = zero - nvelsav = zero - ndfree = zero - nfixat = zero - tdelta = 0.0 - usebox = zero - if (use_bounds) usebox = one - use4d = zero - usefq = zero - merged = zero - vcharmm = 24 - ntitle = one - write (idcd) header,nframe,nprev,ncrdsav,nstep, - & nvelsav,zero,zero,ndfree,nfixat, - & tdelta,usebox,use4d,usefq,merged, - & zero,zero,zero,zero,zero,vcharmm - write (idcd) ntitle,title(1:80) - write (idcd) n - end if -c -c append the lattice values based on header flag value; -c using angle values is NAMD style, cosine values is CHARMM -c - if (use_bounds) then -c write (idcd) xbox,gamma_cos,ybox,beta_cos,alpha_cos,zbox - write (idcd) xbox,gamma,ybox,beta,alpha,zbox - end if -c -c append the force components along each axis in turn -c - write (idcd) (real(-desum(1,i)),i=1,n) - write (idcd) (real(-desum(2,i)),i=1,n) - write (idcd) (real(-desum(3,i)),i=1,n) -c -c close the output unit if opened by this routine -c - if (.not. opened) close (unit=idcd) - return - end diff --git a/source/prtuind.f b/source/prtuind.f deleted file mode 100644 index 485963007..000000000 --- a/source/prtuind.f +++ /dev/null @@ -1,258 +0,0 @@ -c -c -c ################################################### -c ## COPYRIGHT (C) 2023 by Jay William Ponder ## -c ## All Rights Reserved ## -c ################################################### -c -c ################################################################ -c ## ## -c ## subroutine prtuind -- output of atomic induced dipoles ## -c ## ## -c ################################################################ -c -c -c "prtuind" writes out a set of induced dipole components -c to an external file in Tinker XYZ format -c -c - subroutine prtuind (iind) - use atomid - use atoms - use bound - use boxes - use couple - use files - use inform - use polar - use potent - use solpot - use titles - use units - implicit none - integer i,j,k,iind - integer size,crdsiz - real*8 crdmin,crdmax - logical opened - character*2 atmc - character*2 crdc - character*2 digc - character*25 fstr - character*240 indfile -c -c -c open the output unit if not already done -c - inquire (unit=iind,opened=opened) - if (.not. opened) then - indfile = filename(1:leng)//'.uind' - call version (indfile,'new') - open (unit=iind,file=indfile,status='new') - end if -c -c check for large systems needing extended formatting -c - atmc = 'i6' - if (n .ge. 100000) atmc = 'i7' - if (n .ge. 1000000) atmc = 'i8' - crdmin = 0.0d0 - crdmax = 0.0d0 - do i = 1, n - crdmin = min(crdmin,x(i),y(i),z(i)) - crdmax = max(crdmax,x(i),y(i),z(i)) - end do - crdsiz = 6 - if (crdmin .le. -1000.0d0) crdsiz = 7 - if (crdmax .ge. 10000.0d0) crdsiz = 7 - if (crdmin .le. -10000.0d0) crdsiz = 8 - if (crdmax .ge. 100000.0d0) crdsiz = 8 - crdsiz = crdsiz + max(6,digits) - size = 0 - call numeral (crdsiz,crdc,size) - if (digits .le. 6) then - digc = '6 ' - else if (digits .le. 8) then - digc = '8' - else - digc = '10' - end if -c -c write out the number of atoms and the title -c - if (ltitle .eq. 0) then - fstr = '('//atmc//')' - write (iind,fstr(1:4)) n - else - fstr = '('//atmc//',2x,a)' - write (iind,fstr(1:9)) n,title(1:ltitle) - end if -c -c write out the periodic cell lengths and angles -c - if (use_bounds) then - fstr = '(1x,6f'//crdc//'.'//digc//')' - write (iind,fstr) xbox,ybox,zbox,alpha,beta,gamma - end if -c -c write out the induced dipole components for each atom -c - fstr = '('//atmc//',2x,a3,3f'//crdc// - & '.'//digc//',i6,8'//atmc//')' - if (use_solv .and. - & (solvtyp(1:2).eq.'GK'.or.solvtyp(1:2).eq.'PB')) then - do i = 1, n - k = n12(i) - if (k .eq. 0) then - write (iind,fstr) i,name(i),(debye*uinds(j,i),j=1,3), - & type(i) - else - write (iind,fstr) i,name(i),(debye*uinds(j,i),j=1,3), - & type(i),(i12(j,i),j=1,k) - end if - end do - else - do i = 1, n - k = n12(i) - if (k .eq. 0) then - write (iind,fstr) i,name(i),(debye*uind(j,i),j=1,3), - & type(i) - else - write (iind,fstr) i,name(i),(debye*uind(j,i),j=1,3), - & type(i),(i12(j,i),j=1,k) - end if - end do - end if -c -c close the output unit if opened by this routine -c - if (.not. opened) close (unit=iind) - return - end -c -c -c ############################################################# -c ## ## -c ## subroutine prtdcdu -- output of DCD induced dipoles ## -c ## ## -c ############################################################# -c -c -c "prtdcdu" writes out a set of induced dipole components to -c a file in CHARMM DCD binary format compatible with the VMD -c visualization software and other packages -c -c note the format used is based on the "dcdplugin.c" code from -c the NAMD and VMD programs, and tutorial 4.1 from the software -c package GENESIS: Generalized-Ensemble Simulation System -c -c variables and parameters: -c -c header type of data (CORD=coordinates, VELD=velocities) -c nframe number of frames stored in the DCD file -c nprev number of previous integration steps -c ncrdsav frequency in steps for saving coordinate frames -c nstep number of integration steps in the total run -c nvelsav frequency of coordinate saves with velocity data -c ndfree number of degrees of freedom for the system -c nfixat number of fixed atoms for the system -c usebox flag for periodic boundaries (1=true, 0=false) -c use4d flag for 4D trajectory (1=true, 0=false) -c usefq flag for fluctuating charges (1=true, 0=false) -c merged result of merge without checks (1=true, 0=false) -c vcharmm version of CHARMM software for compatibility -c -c in general a value of zero for any of the above indicates that -c the particular feature is unused -c -c - subroutine prtdcdu (idcd,first) - use atoms - use bound - use boxes - use files - use polar - use potent - use solpot - use titles - use units - implicit none - integer i,idcd - integer zero,one - integer nframe,nprev - integer ncrdsav,nstep - integer nvelsav,ndfree - integer nfixat,usebox - integer use4d,usefq - integer merged,vcharmm - integer ntitle - real*4 tdelta - logical opened,first - character*4 header - character*240 dcdfile -c -c -c open the output unit if not already done -c - inquire (unit=idcd,opened=opened) - if (.not. opened) then - dcdfile = filename(1:leng)//'.dcdu' - call version (dcdfile,'new') - open (unit=idcd,file=dcdfile,form='unformatted',status='new') - end if -c -c write header info along with title and number of atoms -c - if (first) then - first = .false. - zero = 0 - one = 1 - header = 'CORD' - nframe = zero - nprev = zero - ncrdsav = one - nstep = zero - nvelsav = zero - ndfree = zero - nfixat = zero - tdelta = 0.0 - usebox = zero - if (use_bounds) usebox = one - use4d = zero - usefq = zero - merged = zero - vcharmm = 24 - ntitle = one - write (idcd) header,nframe,nprev,ncrdsav,nstep, - & nvelsav,zero,zero,ndfree,nfixat, - & tdelta,usebox,use4d,usefq,merged, - & zero,zero,zero,zero,zero,vcharmm - write (idcd) ntitle,title(1:80) - write (idcd) n - end if -c -c append the lattice values based on header flag value; -c using angle values is NAMD style, cosine values is CHARMM -c - if (use_bounds) then -c write (idcd) xbox,gamma_cos,ybox,beta_cos,alpha_cos,zbox - write (idcd) xbox,gamma,ybox,beta,alpha,zbox - end if -c -c append the induced dipoles along each axis in turn -c - if (use_solv .and. - & (solvtyp(1:2).eq.'GK'.or.solvtyp(1:2).eq.'PB')) then - write (idcd) (real(debye*uinds(1,i)),i=1,n) - write (idcd) (real(debye*uinds(2,i)),i=1,n) - write (idcd) (real(debye*uinds(3,i)),i=1,n) - else - write (idcd) (real(debye*uind(1,i)),i=1,n) - write (idcd) (real(debye*uind(2,i)),i=1,n) - write (idcd) (real(debye*uind(3,i)),i=1,n) - end if -c -c close the output unit if opened by this routine -c - if (.not. opened) close (unit=idcd) - return - end diff --git a/source/prtvel.f b/source/prtvel.f deleted file mode 100644 index 0b9608cbb..000000000 --- a/source/prtvel.f +++ /dev/null @@ -1,230 +0,0 @@ -c -c -c ################################################### -c ## COPYRIGHT (C) 2023 by Jay William Ponder ## -c ## All Rights Reserved ## -c ################################################### -c -c ############################################################ -c ## ## -c ## subroutine prtvel -- output of velocity components ## -c ## ## -c ############################################################ -c -c -c "prtvel" writes out a set of atomic velocity components -c to an external file in Tinker XYZ format -c -c - subroutine prtvel (ivel) - use atomid - use atoms - use bound - use boxes - use couple - use files - use inform - use moldyn - use titles - implicit none - integer i,j,k,ivel - integer size,crdsiz - real*8 crdmin,crdmax - logical opened - character*2 atmc - character*2 crdc - character*2 digc - character*25 fstr - character*240 velfile -c -c -c open the output unit if not already done -c - inquire (unit=ivel,opened=opened) - if (.not. opened) then - velfile = filename(1:leng)//'.vel' - call version (velfile,'new') - open (unit=ivel,file=velfile,status='new') - end if -c -c check for large systems needing extended formatting -c - atmc = 'i6' - if (n .ge. 100000) atmc = 'i7' - if (n .ge. 1000000) atmc = 'i8' - crdmin = 0.0d0 - crdmax = 0.0d0 - do i = 1, n - crdmin = min(crdmin,x(i),y(i),z(i)) - crdmax = max(crdmax,x(i),y(i),z(i)) - end do - crdsiz = 6 - if (crdmin .le. -1000.0d0) crdsiz = 7 - if (crdmax .ge. 10000.0d0) crdsiz = 7 - if (crdmin .le. -10000.0d0) crdsiz = 8 - if (crdmax .ge. 100000.0d0) crdsiz = 8 - crdsiz = crdsiz + max(6,digits) - size = 0 - call numeral (crdsiz,crdc,size) - if (digits .le. 6) then - digc = '6 ' - else if (digits .le. 8) then - digc = '8' - else - digc = '10' - end if -c -c write out the number of atoms and the title -c - if (ltitle .eq. 0) then - fstr = '('//atmc//')' - write (ivel,fstr(1:4)) n - else - fstr = '('//atmc//',2x,a)' - write (ivel,fstr(1:9)) n,title(1:ltitle) - end if -c -c write out the periodic cell lengths and angles -c - if (use_bounds) then - fstr = '(1x,6f'//crdc//'.'//digc//')' - write (ivel,fstr) xbox,ybox,zbox,alpha,beta,gamma - end if -c -c write out the velocity components for each atom -c - fstr = '('//atmc//',2x,a3,3f'//crdc// - & '.'//digc//',i6,8'//atmc//')' - do i = 1, n - k = n12(i) - if (k .eq. 0) then - write (ivel,fstr) i,name(i),(v(j,i),j=1,3),type(i) - else - write (ivel,fstr) i,name(i),(v(j,i),j=1,3),type(i), - & (i12(j,i),j=1,k) - end if - end do -c -c close the output unit if opened by this routine -c - if (.not. opened) close (unit=ivel) - return - end -c -c -c ################################################################# -c ## ## -c ## subroutine prtdcdv -- output of DCD velocity components ## -c ## ## -c ################################################################# -c -c -c "prtdcdv" writes out a set of atomic velocity components to -c a file in CHARMM DCD binary format compatible with the VMD -c visualization software and other packages -c -c note the format used is based on the "dcdplugin.c" code from -c the NAMD and VMD programs, and tutorial 4.1 from the software -c package GENESIS: Generalized-Ensemble Simulation System -c -c variables and parameters: -c -c header type of data (CORD=coordinates, VELD=velocities) -c nframe number of frames stored in the DCD file -c nprev number of previous integration steps -c ncrdsav frequency in steps for saving coordinate frames -c nstep number of integration steps in the total run -c nvelsav frequency of coordinate saves with velocity data -c ndfree number of degrees of freedom for the system -c nfixat number of fixed atoms for the system -c usebox flag for periodic boundaries (1=true, 0=false) -c use4d flag for 4D trajectory (1=true, 0=false) -c usefq flag for fluctuating charges (1=true, 0=false) -c merged result of merge without checks (1=true, 0=false) -c vcharmm version of CHARMM software for compatibility -c -c in general a value of zero for any of the above indicates that -c the particular feature is unused -c -c - subroutine prtdcdv (idcd,first) - use atoms - use bound - use boxes - use files - use moldyn - use titles - implicit none - integer i,idcd - integer zero,one - integer nframe,nprev - integer ncrdsav,nstep - integer nvelsav,ndfree - integer nfixat,usebox - integer use4d,usefq - integer merged,vcharmm - integer ntitle - real*4 tdelta - logical opened,first - character*4 header - character*240 dcdfile -c -c -c open the output unit if not already done -c - inquire (unit=idcd,opened=opened) - if (.not. opened) then - dcdfile = filename(1:leng)//'.dcdv' - call version (dcdfile,'new') - open (unit=idcd,file=dcdfile,form='unformatted',status='new') - end if -c -c write header info along with title and number of atoms -c - if (first) then - first = .false. - zero = 0 - one = 1 - header = 'CORD' - nframe = zero - nprev = zero - ncrdsav = one - nstep = zero - nvelsav = zero - ndfree = zero - nfixat = zero - tdelta = 0.0 - usebox = zero - if (use_bounds) usebox = one - use4d = zero - usefq = zero - merged = zero - vcharmm = 24 - ntitle = one - write (idcd) header,nframe,nprev,ncrdsav,nstep, - & nvelsav,zero,zero,ndfree,nfixat, - & tdelta,usebox,use4d,usefq,merged, - & zero,zero,zero,zero,zero,vcharmm - write (idcd) ntitle,title(1:80) - write (idcd) n - end if -c -c append the lattice values based on header flag value; -c using angle values is NAMD style, cosine values is CHARMM -c - if (use_bounds) then -c write (idcd) xbox,gamma_cos,ybox,beta_cos,alpha_cos,zbox - write (idcd) xbox,gamma,ybox,beta,alpha,zbox - end if -c -c append the velocity components along each axis in turn -c - write (idcd) (real(v(1,i)),i=1,n) - write (idcd) (real(v(2,i)),i=1,n) - write (idcd) (real(v(3,i)),i=1,n) -c -c close the output unit if opened by this routine -c - if (.not. opened) close (unit=idcd) - return - end diff --git a/source/prtxyz.f b/source/prtxyz.f index e1b19106d..1a3101e6d 100644 --- a/source/prtxyz.f +++ b/source/prtxyz.f @@ -24,9 +24,11 @@ subroutine prtxyz (ixyz) use couple use files use inform + use output use titles implicit none integer i,j,k,ixyz + integer ii integer size,crdsiz real*8 crdmin,crdmax logical opened @@ -75,12 +77,22 @@ subroutine prtxyz (ixyz) c c write out the number of atoms and the title c - if (ltitle .eq. 0) then - fstr = '('//atmc//')' - write (ixyz,fstr(1:4)) n + if (.not. onlysave) then + if (ltitle .eq. 0) then + fstr = '('//atmc//')' + write (ixyz,fstr(1:4)) n + else + fstr = '('//atmc//',2x,a)' + write (ixyz,fstr(1:9)) n,title(1:ltitle) + end if else - fstr = '('//atmc//',2x,a)' - write (ixyz,fstr(1:9)) n,title(1:ltitle) + if (ltitle .eq. 0) then + fstr = '('//atmc//')' + write (ixyz,fstr(1:4)) nonly + else + fstr = '('//atmc//',2x,a)' + write (ixyz,fstr(1:9)) nonly,title(1:ltitle) + end if end if c c write out the periodic cell lengths and angles @@ -94,19 +106,371 @@ subroutine prtxyz (ixyz) c fstr = '('//atmc//',2x,a3,3f'//crdc// & '.'//digc//',i6,8'//atmc//')' + if (.not. onlysave) then + do i = 1, n + k = n12(i) + if (k .eq. 0) then + write (ixyz,fstr) i,name(i),x(i),y(i),z(i),type(i) + else + write (ixyz,fstr) i,name(i),x(i),y(i),z(i),type(i), + & (i12(j,i),j=1,k) + end if + end do + else + do ii = 1, nonly + i = ionly(ii) + k = n12(i) + if (k .eq. 0) then + write (ixyz,fstr) ii,name(i),x(i),y(i),z(i),type(i) + else + write (ixyz,fstr) ii,name(i),x(i),y(i),z(i),type(i), + & (ionlyinv(i12(j,i)),j=1,k) + end if + end do + end if +c +c close the output unit if opened by this routine +c + if (.not. opened) close (unit=ixyz) + return + end +c +c +c ########################################################### +c ## ## +c ## subroutine prtvec3 -- output of atomic 3D vectors ## +c ## ## +c ########################################################### +c +c +c "prtvec3" writes out a set of atomic vectors +c to an external disk file in Tinker XYZ format +c +c + subroutine prtvec3 (iunit,mode) + use atomid + use atoms + use bound + use boxes + use couple + use files + use inform + use output + use titles + implicit none + integer i,j,k + integer ii,iunit + integer size,crdsiz + real*8 crdmin,crdmax + logical opened + character*2 atmc + character*2 crdc + character*2 digc + character*3 mode + character*25 fstr + character*240 outputfile +c +c +c open the output unit if not already done +c + inquire (unit=iunit,opened=opened) + if (.not. opened) then + if (mode .eq. 'XYZ') then + outputfile = filename(1:leng)//'.xyz' + else if (mode .eq. 'VEL') then + outputfile = filename(1:leng)//'.vel' + else if (mode .eq. 'FRC') then + outputfile = filename(1:leng)//'.frc' + else if (mode .eq. 'UIN') then + outputfile = filename(1:leng)//'.uind' + else if (mode .eq. 'UST') then + outputfile = filename(1:leng)//'.ustc' + else if (mode .eq. 'UCH') then + outputfile = filename(1:leng)//'.uchg' + else if (mode .eq. 'UDR') then + outputfile = filename(1:leng)//'.udir' + else if (mode .eq. 'DEF') then + outputfile = filename(1:leng)//'.def' + else if (mode .eq. 'TEF') then + outputfile = filename(1:leng)//'.tef' + end if + call version (outputfile,'new') + open (unit=iunit,file=outputfile,status='new') + end if +c +c copy the vector values to a common array +c + call copyvec3 (print3n,mode) +c +c check for large systems needing extended formatting +c + atmc = 'i6' + if (n .ge. 100000) atmc = 'i7' + if (n .ge. 1000000) atmc = 'i8' + crdmin = 0.0d0 + crdmax = 0.0d0 do i = 1, n - k = n12(i) - if (k .eq. 0) then - write (ixyz,fstr) i,name(i),x(i),y(i),z(i),type(i) + crdmin = min(crdmin,print3n(1,i),print3n(2,i),print3n(3,i)) + crdmax = max(crdmax,print3n(1,i),print3n(2,i),print3n(3,i)) + end do + crdsiz = 6 + if (crdmin .le. -1000.0d0) crdsiz = 7 + if (crdmax .ge. 10000.0d0) crdsiz = 7 + if (crdmin .le. -10000.0d0) crdsiz = 8 + if (crdmax .ge. 100000.0d0) crdsiz = 8 + crdsiz = crdsiz + max(6,digits) + size = 0 + call numeral (crdsiz,crdc,size) + if (digits .le. 6) then + digc = '6 ' + else if (digits .le. 8) then + digc = '8' + else + digc = '10' + end if +c +c write out the number of atoms and the title +c + if (.not. onlysave) then + if (ltitle .eq. 0) then + fstr = '('//atmc//')' + write (iunit,fstr(1:4)) n else - write (ixyz,fstr) i,name(i),x(i),y(i),z(i),type(i), - & (i12(j,i),j=1,k) + fstr = '('//atmc//',2x,a)' + write (iunit,fstr(1:9)) n,title(1:ltitle) end if - end do + else + if (ltitle .eq. 0) then + fstr = '('//atmc//')' + write (iunit,fstr(1:4)) nonly + else + fstr = '('//atmc//',2x,a)' + write (iunit,fstr(1:9)) nonly,title(1:ltitle) + end if + end if +c +c write out the periodic cell lengths and angles +c + if (use_bounds) then + fstr = '(1x,6f'//crdc//'.'//digc//')' + write (iunit,fstr) xbox,ybox,zbox,alpha,beta,gamma + end if +c +c write out the atomic coordinates for each atom +c + fstr = '('//atmc//',2x,a3,3f'//crdc// + & '.'//digc//',i6,8'//atmc//')' + if (.not. onlysave) then + do i = 1, n + k = n12(i) + if (k .eq. 0) then + write (iunit,fstr) i,name(i),(print3n(j,i),j=1,3), + & type(i) + else + write (iunit,fstr) i,name(i),(print3n(j,i),j=1,3), + & type(i),(i12(j,i),j=1,k) + end if + end do + else + do ii = 1, nonly + i = ionly(ii) + k = n12(i) + if (k .eq. 0) then + write (iunit,fstr) ii,name(i),(print3n(j,i),j=1,3), + & type(i) + else + write (iunit,fstr) ii,name(i),(print3n(j,i),j=1,3), + & type(i),(ionlyinv(i12(j,i)),j=1,k) + end if + end do + end if c c close the output unit if opened by this routine c - if (.not. opened) close (unit=ixyz) + if (.not. opened) close (unit=iunit) + return + end +c +c +c ############################################################## +c ## ## +c ## subroutine copyvec3 -- copy vector values to print3n ## +c ## ## +c ############################################################## +c +c +c "copyvec3" copies the values of the requested vector to the +c common array "print3n" +c +c + subroutine copyvec3 (print3n,mode) + use atoms + use boxes + use dipole + use deriv + use charge + use expol + use moldyn + use mpole + use polar + use polpot + use potent + use solpot + use units + implicit none + integer i,j,k + integer ii + real*8 c + real*8 xi,yi,zi,ri + real*8 print3n(3,*) + character*3 mode +c +c +c copy the requested vector values to the common array "print3n" +c + if (mode .eq. 'XYZ') then + do i = 1, n + print3n(1,i) = x(i) + print3n(2,i) = y(i) + print3n(3,i) = z(i) + end do + else if (mode .eq. 'VEL') then + do i = 1, n + print3n(1,i) = v(1,i) + print3n(2,i) = v(2,i) + print3n(3,i) = v(3,i) + end do + else if (mode .eq. 'FRC') then + do i = 1, n + print3n(1,i) = -desum(1,i) + print3n(2,i) = -desum(2,i) + print3n(3,i) = -desum(3,i) + end do + else if (mode .eq. 'UIN') then + if (use_solv .and. + & (solvtyp(1:2).eq.'GK'.or.solvtyp(1:2).eq.'PB')) then + do i = 1, n + print3n(1,i) = debye*uinds(1,i) + print3n(2,i) = debye*uinds(2,i) + print3n(3,i) = debye*uinds(3,i) + end do + else + do i = 1, n + print3n(1,i) = debye*uind(1,i) + print3n(2,i) = debye*uind(2,i) + print3n(3,i) = debye*uind(3,i) + end do + end if + else if (mode .eq. 'UST') then + do i = 1, ndipole + j = idpl(1,i) + k = idpl(2,i) + xi = x(j) - x(k) + yi = y(j) - y(k) + zi = z(j) - z(k) + ri = sqrt(xi*xi + yi*yi + zi*zi) + print3n(1,i) = bdpl(i) * (xi/ri) + print3n(2,i) = bdpl(i) * (yi/ri) + print3n(3,i) = bdpl(i) * (zi/ri) + end do + do ii = 1, npole + i = ipole(ii) + print3n(1,i) = debye*rpole(2,i) + print3n(2,i) = debye*rpole(3,i) + print3n(3,i) = debye*rpole(4,i) + end do + else if (mode .eq. 'UCH') then + do ii = 1, nion + i = iion(ii) + c = pchg(i) + print3n(1,i) = c*debye*(x(i)-xcenter) + print3n(2,i) = c*debye*(y(i)-ycenter) + print3n(3,i) = c*debye*(z(i)-zcenter) + end do + do ii = 1, npole + i = ipole(ii) + c = rpole(1,i) + print3n(1,i) = c*debye*(x(i)-xcenter) + print3n(2,i) = c*debye*(y(i)-ycenter) + print3n(3,i) = c*debye*(z(i)-zcenter) + end do + else if (mode .eq. 'UDR') then + if (.not. use_expol) then + if (use_solv .and. + & (solvtyp(1:2).eq.'GK'.or.solvtyp(1:2).eq.'PB')) then + do i = 1, n + do j = 1, 3 + print3n(j,i) = debye*udirs(j,i) + end do + end do + else + do i = 1, n + do j = 1, 3 + print3n(j,i) = debye*udir(j,i) + end do + end do + end if + else + do i = 1, n + do j = 1, 3 + print3n(j,i) = 0.0d0 + do k = 1, 3 + print3n(j,i) = print3n(j,i) + & + udir(k,i)*polinv(j,k,i) + end do + print3n(j,i) = debye*print3n(j,i) + end do + end do + end if + else if (mode .eq. 'DEF') then + if (use_solv .and. + & (solvtyp(1:2).eq.'GK'.or.solvtyp(1:2).eq.'PB')) then + do i = 1, n + c = elefield/polarity(i) + do j = 1, 3 + print3n(j,i) = c * udirs(j,i) + end do + end do + else + do i = 1, n + c = elefield/polarity(i) + do j = 1, 3 + print3n(j,i) = c * udir(j,i) + end do + end do + end if + else if (mode .eq. 'TEF') then + if (.not. use_expol) then + if (use_solv .and. + & (solvtyp(1:2).eq.'GK'.or.solvtyp(1:2).eq.'PB')) then + do i = 1, n + c = elefield/polarity(i) + do j = 1, 3 + print3n(j,i) = c * uinds(j,i) + end do + end do + else + do i = 1, n + c = elefield/polarity(i) + do j = 1, 3 + print3n(j,i) = c * uind(j,i) + end do + end do + end if + else + do i = 1, n + c = elefield/polarity(i) + do j = 1, 3 + print3n(j,i) = 0.0d0 + do k = 1, 3 + print3n(j,i) = print3n(j,i) + & + uind(k,i)*polscale(j,k,i) + end do + print3n(j,i) = c * print3n(j,i) + end do + end do + end if + end if return end c @@ -151,6 +515,7 @@ subroutine prtdcd (idcd,first) use bound use boxes use files + use output use titles implicit none integer i,idcd @@ -204,7 +569,11 @@ subroutine prtdcd (idcd,first) & tdelta,usebox,use4d,usefq,merged, & zero,zero,zero,zero,zero,vcharmm write (idcd) ntitle,title(1:80) - write (idcd) n + if (.not. onlysave) then + write (idcd) n + else + write (idcd) nonly + end if end if c c append the lattice values based on header flag value; @@ -217,12 +586,174 @@ subroutine prtdcd (idcd,first) c c append the atomic coordinates along each axis in turn c - write (idcd) (real(x(i)),i=1,n) - write (idcd) (real(y(i)),i=1,n) - write (idcd) (real(z(i)),i=1,n) + if (.not. onlysave) then + write (idcd) (real(x(i)),i=1,n) + write (idcd) (real(y(i)),i=1,n) + write (idcd) (real(z(i)),i=1,n) + else + write (idcd) (real(x(ionly(i))),i=1,nonly) + write (idcd) (real(y(ionly(i))),i=1,nonly) + write (idcd) (real(z(ionly(i))),i=1,nonly) + end if c c close the output unit if opened by this routine c if (.not. opened) close (unit=idcd) return end +c +c +c ################################################################ +c ## ## +c ## subroutine prtdcdv3 -- output of DCD atomic 3D vectors ## +c ## ## +c ################################################################ +c +c +c "prtdcdV3" writes out a set of atomic 3D vectors to +c a file in CHARMM DCD binary format compatible with the VMD +c visualization software and other packages +c +c note the format used is based on the "dcdplugin.c" code from +c the NAMD and VMD programs, and tutorial 4.1 from the software +c package GENESIS: Generalized-Ensemble Simulation System +c +c variables and parameters: +c +c header type of data (CORD=coordinates, VELD=velocities) +c nframe number of frames stored in the DCD file +c nprev number of previous integration steps +c ncrdsav frequency in steps for saving coordinate frames +c nstep number of integration steps in the total run +c nvelsav frequency of coordinate saves with velocity data +c ndfree number of degrees of freedom for the system +c nfixat number of fixed atoms for the system +c usebox flag for periodic boundaries (1=true, 0=false) +c use4d flag for 4D trajectory (1=true, 0=false) +c usefq flag for fluctuating charges (1=true, 0=false) +c merged result of merge without checks (1=true, 0=false) +c vcharmm version of CHARMM software for compatibility +c +c in general a value of zero for any of the above indicates that +c the particular feature is unused +c +c + subroutine prtdcdv3 (idcd,first,mode) + use atoms + use bound + use boxes + use files + use output + use titles + implicit none + integer i,idcd + integer zero,one + integer nframe,nprev + integer ncrdsav,nstep + integer nvelsav,ndfree + integer nfixat,usebox + integer use4d,usefq + integer merged,vcharmm + integer ntitle + real*4 tdelta + logical opened,first + character*3 mode + character*4 header + character*240 dcdfile +c +c +c open the output unit if not already done +c + inquire (unit=idcd,opened=opened) + if (.not. opened) then + dcdfile = filename(1:leng)//'.dcd' + call version (dcdfile,'new') + open (unit=idcd,file=dcdfile,form='unformatted',status='new') + end if + if (.not. opened) then + if (mode .eq. 'XYZ') then + dcdfile = filename(1:leng)//'.dcd' + else if (mode .eq. 'VEL') then + dcdfile = filename(1:leng)//'.dcdv' + else if (mode .eq. 'FRC') then + dcdfile = filename(1:leng)//'.dcdf' + else if (mode .eq. 'UIN') then + dcdfile = filename(1:leng)//'.dcdui' + else if (mode .eq. 'UST') then + dcdfile = filename(1:leng)//'.dcdus' + else if (mode .eq. 'UCH') then + dcdfile = filename(1:leng)//'.dcduc' + else if (mode .eq. 'UDR') then + dcdfile = filename(1:leng)//'.dcdud' + else if (mode .eq. 'DEF') then + dcdfile = filename(1:leng)//'.dcdde' + else if (mode .eq. 'TEF') then + dcdfile = filename(1:leng)//'.dcdte' + end if + call version (dcdfile,'new') + open (unit=idcd,file=dcdfile,form='unformatted',status='new') + end if +c +c write header info along with title and number of atoms +c + if (first) then + first = .false. + zero = 0 + one = 1 + header = 'CORD' + nframe = zero + nprev = zero + ncrdsav = one + nstep = zero + nvelsav = zero + ndfree = zero + nfixat = zero + tdelta = 0.0 + usebox = zero + if (use_bounds) usebox = one + use4d = zero + usefq = zero + merged = zero + vcharmm = 24 + ntitle = one + write (idcd) header,nframe,nprev,ncrdsav,nstep, + & nvelsav,zero,zero,ndfree,nfixat, + & tdelta,usebox,use4d,usefq,merged, + & zero,zero,zero,zero,zero,vcharmm + write (idcd) ntitle,title(1:80) + if (.not. onlysave) then + write (idcd) n + else + write (idcd) nonly + end if + end if +c +c append the lattice values based on header flag value; +c using angle values is NAMD style, cosine values is CHARMM +c + if (use_bounds) then +c write (idcd) xbox,gamma_cos,ybox,beta_cos,alpha_cos,zbox + write (idcd) xbox,gamma,ybox,beta,alpha,zbox + end if +c +c copy the vector values to a common array +c + call copyvec3 (print3n,mode) +c +c append the atomic coordinates along each axis in turn +c + if (.not. onlysave) then + write (idcd) (real(print3n(1,i)),i=1,n) + write (idcd) (real(print3n(2,i)),i=1,n) + write (idcd) (real(print3n(3,i)),i=1,n) + else + write (idcd) (real(print3n(1,ionly(i))),i=1,nonly) + write (idcd) (real(print3n(2,ionly(i))),i=1,nonly) + write (idcd) (real(print3n(3,ionly(i))),i=1,nonly) + end if +c +c close the output unit if opened by this routine +c + if (.not. opened) close (unit=idcd) + return + end \ No newline at end of file diff --git a/source/uatom.f b/source/uatom.f new file mode 100644 index 000000000..ea7c74eb2 --- /dev/null +++ b/source/uatom.f @@ -0,0 +1,33 @@ +c +c +c ############################################################## +c ## COPYRIGHT (C) 2024 by Moses KJ Chung & Jay W Ponder ## +c ## All Rights Reserved ## +c ############################################################## +c +c ####################################################### +c ## ## +c ## module uatom -- properties based on atom type ## +c ## ## +c ####################################################### +c +c +c nunique number of unique atom types in the system +c utype map from unique type to atom type +c utypeinv map from atom type to unique type +c utv1 unique type vector 1 +c utv2 unique type vector 2 +c utv3 unique type vector 3 +c +c + module uatom + use sizes + implicit none + integer nunique + integer utype(maxtyp) + integer utypeinv(maxtyp) + real*8 utv1(3,maxtyp) + real*8 utv2(3,maxtyp) + real*8 utv3(3,maxtyp) + save + end diff --git a/source/uniquetyp.f b/source/uniquetyp.f new file mode 100644 index 000000000..0e3b2c868 --- /dev/null +++ b/source/uniquetyp.f @@ -0,0 +1,101 @@ +c +c +c ############################################################## +c ## COPYRIGHT (C) 2024 by Moses KJ Chung & Jay W Ponder ## +c ## All Rights Reserved ## +c ############################################################## +c +c ################################################################## +c ## ## +c ## subroutine uniquetyp -- get unique atom type information ## +c ## ## +c ################################################################## +c +c +c "uniquetyp" determines the number of unique types and map between +c atom types and unique types +c +c + subroutine uniquetyp + use atoms + use uatom + implicit none + integer i,j + integer at + real*8 xx +c +c +c initialize nunique, utype, and utypeinv +c + nunique = 0 + do i = 1, maxtyp + utype(i) = 0 + utypeinv(i) = 0 + end do +c +c loop through atoms to find unique atom types +c + do i = 1, n + at = type(i) +c +c check if the current atom type is already in unique type array +c + do j = 1, nunique + if (at .eq. utype(j)) goto 10 + end do +c +c set unique atom types +c + nunique = nunique + 1 + utype(nunique) = at + utypeinv(at) = nunique +10 continue + end do + return + end +c +c +c ########################################################### +c ## ## +c ## subroutine velunique -- unique atom type velocity ## +c ## ## +c ########################################################### +c +c +c "velunique" computes the unique atom type velocities +c +c + subroutine velunique + use atomid + use atoms + use moldyn + use output + use uatom + implicit none + integer i + integer ut +c +c +c zero out velocity +c + do i = 1, nunique + utv1(1,i) = 0.0d0 + utv1(2,i) = 0.0d0 + utv1(3,i) = 0.0d0 + end do +!$OMP PARALLEL default(private) +!$OMP& shared(n,utv1,type,utypeinv,v) +!$OMP DO reduction(+:utv1) +c +c compute velocity by unique atom type +c + do i = 1, n + ut = utypeinv(type(i)) + utv1(1,ut) = utv1(1,ut) + v(1,i) + utv1(2,ut) = utv1(2,ut) + v(2,i) + utv1(3,ut) = utv1(3,ut) + v(3,i) + end do +!$OMP END DO +!$OMP END PARALLEL + return + end diff --git a/windows/cygwin/compile.make b/windows/cygwin/compile.make index 564a3c3e5..32a8a4a71 100755 --- a/windows/cygwin/compile.make +++ b/windows/cygwin/compile.make @@ -528,14 +528,11 @@ gfortran -c -O3 -fopenmp prtarc.f gfortran -c -O3 -fopenmp prtcif.f gfortran -c -O3 -fopenmp prtdyn.f gfortran -c -O3 -fopenmp prterr.f -gfortran -c -O3 -fopenmp prtfrc.f gfortran -c -O3 -fopenmp prtint.f gfortran -c -O3 -fopenmp prtmol2.f gfortran -c -O3 -fopenmp prtpdb.f gfortran -c -O3 -fopenmp prtprm.f gfortran -c -O3 -fopenmp prtseq.f -gfortran -c -O3 -fopenmp prtuind.f -gfortran -c -O3 -fopenmp prtvel.f gfortran -c -O3 -fopenmp prtxyz.f gfortran -c -O3 -fopenmp pss.f gfortran -c -O3 -fopenmp pssrigid.f diff --git a/windows/cygwin/library.make b/windows/cygwin/library.make index 33f750b0b..b6757e4ad 100755 --- a/windows/cygwin/library.make +++ b/windows/cygwin/library.make @@ -438,14 +438,11 @@ prtarc.o \ prtcif.o \ prtdyn.o \ prterr.o \ -prtfrc.o \ prtint.o \ prtmol2.o \ prtpdb.o \ prtprm.o \ prtseq.o \ -prtuind.o \ -prtvel.o \ prtxyz.o \ ptable.o \ qmstuf.o \ diff --git a/windows/intel/compile.bat b/windows/intel/compile.bat index 62788ac1e..d6b24840a 100755 --- a/windows/intel/compile.bat +++ b/windows/intel/compile.bat @@ -529,14 +529,11 @@ ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtarc.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtcif.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtdyn.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prterr.f -ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtfrc.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtint.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtmol2.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtpdb.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtprm.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtseq.f -ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtuind.f -ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtvel.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp prtxyz.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp pss.f ifort /c /O3 /QxHost /Qip- /Qprec-div- /w /Qopenmp pssrigid.f diff --git a/windows/intel/generic.bat b/windows/intel/generic.bat index e79a0264d..6ca9c5a5d 100755 --- a/windows/intel/generic.bat +++ b/windows/intel/generic.bat @@ -527,14 +527,11 @@ ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtarc.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtcif.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtdyn.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prterr.f -ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtfrc.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtint.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtmol2.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtpdb.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtprm.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtseq.f -ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtuind.f -ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtvel.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp prtxyz.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp pss.f ifort /c /O3 /arch:sse3 /Qip- /Qprec-div- /w /Qopenmp pssrigid.f diff --git a/windows/intel/library.lbc b/windows/intel/library.lbc index 27480ccd3..51fffd6ce 100644 --- a/windows/intel/library.lbc +++ b/windows/intel/library.lbc @@ -428,14 +428,11 @@ prtarc.obj prtcif.obj prtdyn.obj prterr.obj -prtfrc.obj prtint.obj prtmol2.obj prtpdb.obj prtprm.obj prtseq.obj -prtuind.obj -prtvel.obj prtxyz.obj ptable.obj qmstuf.obj