Removed spurious file changes not belonging to this branch.

variable number	equal 2

variable ofreq	equal 1000

variable efreq	equal 1000

units lj

dimension 3

newton off

boundary  p p p

atom_style hybrid bond ellipsoid

atom_modify sort 0 1.0

# Pair interactions require lists of neighbours to be calculated

neighbor 3.0 bin

neigh_modify every 30 delay 0 check yes

read_data data.duplex2

set atom * mass 3.1575

group all type 1 4

# oxDNA bond interactions - FENE backbone

bond_style oxdna2/fene

bond_coeff * 2.0 0.25 0.7564

# oxDNA pair interactions

pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh

pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32

pair_coeff * * oxdna2/stk     seqav 0.1 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65   

pair_coeff * * oxdna2/hbond   seqav 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45

pair_coeff 1 4 oxdna2/hbond   seqav 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45

pair_coeff 2 3 oxdna2/hbond   seqav 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45

pair_coeff * * oxdna2/xstk    47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 

pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793

pair_coeff * * oxdna2/dh      0.1 1.0 0.815

# NVE ensemble

#fix 1 all nve/dot

fix 1 all nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10

#fix 1 all nve/asphere

#fix 2 all langevin 0.1 0.1 0.03 457145 angmom 10

timestep 1e-5

#comm_style tiled

#fix 3 all balance 10000 1.1 rcb

#compute mol all chunk/atom molecule

#compute mychunk all vcm/chunk mol

#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector

#dump pos all xyz ${ofreq} traj.${number}.xyz

#compute quat all property/atom quatw quati quatj quatk

#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]

#dump_modify quat sort id

#dump_modify quat format line "%d  %13.6le  %13.6le  %13.6le  %13.6le"

compute erot all erotate/asphere

compute ekin all ke

compute epot all pe

variable erot equal c_erot

variable ekin equal c_ekin

variable epot equal c_epot

variable etot equal c_erot+c_ekin+c_epot

fix 5 all print ${efreq} "$(step)  ekin = ${ekin} |  erot = ${erot} | epot = ${epot} | etot = ${etot}" screen no file oxdna_thermo.dat

thermo 100000

dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz

dump_modify out sort id

# dump_modify out format line "%d   %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le  %13.6le %13.6le %13.6le"

run 10000

#write_restart config.${number}.*

/* -*- c++ -*- ----------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.

------------------------------------------------------------------------- */

#ifdef PAIR_CLASS

PairStyle(patchysphere,PairPatchySphere)

#else

#ifndef LMP_PAIR_PATCHY_SPHERE_H

#define LMP_PAIR_PATCHY_SPHERE_H

#include "pair.h"

namespace LAMMPS_NS {

class PairPatchySphere : public Pair {

 public:

  PairPatchySphere(LAMMPS *);

  virtual ~PairPatchySphere();

  virtual void compute(int, int);

  virtual void settings(int, char **);

  void coeff(int, char **);

  virtual void init_style();

  double init_one(int, int);

	//void write_restart(FILE *);

	//void read_restart(FILE *);

	//void write_restart_settings(FILE *);

	//void read_restart_settings(FILE *);

	//void write_data(FILE *);

	//void write_data_all(FILE *);

 protected:

  // enumerate the "encounter" possibilities.

  enum{SPHERE_SPHERE,SPHERE_PATCHY,PATCHY_SPHERE,PATCHY_PATCHY};

  double cut_global;

  double **cut;

  double **inner_cut;

  double **offset;

  double **epsilon, **sigma;

  // Whole bunch of stuff related to bond vectors...

  double **epsilon_b; // Patch binding strength.

  int n_bond_pairs;   // Number of bond pairs (combinations of alpha, beta, ...)

  int **bond_pairs;   // Enumerate bond pairs.

  double **phi_m, **theta_m; // Theta and phi scale factors.

  int n_bonds;            // Number of bond vectors.

  double **bond_vectors;  // Bond vectors.

  class AtomVecEllipsoid *avec;

  bool debug_output;

  void allocate();

  double bond_vec_analytic(const int i, const int j, double a1[3][3],

                           double a2[3][3], double *r12, const double rsq,

                           double *fforce, double *ttor, double *rtor);

  double bond_vec_theta_part(const int i, const int j,

                             const int itype, const int jtype,

                             const double *r12, const double *b_alp,

                             const double *b_bet, double *fforce, double *ttor,

                             double *rtor);

  double bond_vec_phi_part(const int i, const int j,

                           const int itype, const int jtype,

                           const double *r12, const double rsq,

                           const double *b_alp, const double *b_bet,

                           const double *b_gam, const double *b_eps,

                           double *fforce, double *ttor, double *rtor);

};

} // namespace

#endif // LMP_PAIR_PATCHY_SPHERE_H

#endif // PAIR_CLASS

# kokkos_cuda_mpi = KOKKOS/CUDA package, MPICH or OpenMPI with nvcc compiler, Kepler GPU

SHELL = /bin/sh

# ---------------------------------------------------------------------

# compiler/linker settings

# specify flags and libraries needed for your compiler

KOKKOS_ABSOLUTE_PATH = /home/stefan/projects/lammps-mine/lib/kokkos

export MPICH_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper

export OMPI_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper

CC =		mpicxx

CCFLAGS =	-g -O3

SHFLAGS =	-fPIC

DEPFLAGS =	-M

LINK =		mpicxx

LINKFLAGS =	-g -O3

LIB = 

SIZE =		size

ARCHIVE =	ar

ARFLAGS =	-rc

SHLIBFLAGS =	-shared

KOKKOS_DEVICES = Cuda

KOKKOS_ARCH = Kepler35

# ---------------------------------------------------------------------

# LAMMPS-specific settings, all OPTIONAL

# specify settings for LAMMPS features you will use

# if you change any -D setting, do full re-compile after "make clean"

# LAMMPS ifdef settings

# see possible settings in Section 2.2 (step 4) of manual

LMP_INC =	-DLAMMPS_GZIP -DLAMMPS_PNG -DLAMMPS_JPEG -DLAMMPS_FFMPEG

# MPI library

# see discussion in Section 2.2 (step 5) of manual

# MPI wrapper compiler/linker can provide this info

# can point to dummy MPI library in src/STUBS as in Makefile.serial

# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts

# INC = path for mpi.h, MPI compiler settings

# PATH = path for MPI library

# LIB = name of MPI library

MPI_INC =       -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX=1

MPI_PATH = 

MPI_LIB =	

# FFT library

# see discussion in Section 2.2 (step 6) of manaul

# can be left blank to use provided KISS FFT library

# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings

# PATH = path for FFT library

# LIB = name of FFT library

FFT_INC =    	-I/usr/include/

FFT_PATH = 	-L/usr/lib/

FFT_LIB =	-lfftw3

# JPEG and/or PNG library

# see discussion in Section 2.2 (step 7) of manual

# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC

# INC = path(s) for jpeglib.h and/or png.h

# PATH = path(s) for JPEG library and/or PNG library

# LIB = name(s) of JPEG library and/or PNG library

JPG_INC =       

JPG_PATH = 	

JPG_LIB =	-ljpeg -lpng

# ---------------------------------------------------------------------

# build rules and dependencies

# do not edit this section

include	Makefile.package.settings

include	Makefile.package

EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC)

EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH)

EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB)

EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS)

EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS)

# Path to src files

vpath %.cpp ..

vpath %.h ..

# Link target

$(EXE):	$(OBJ) $(EXTRA_LINK_DEPENDS)

	$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)

	$(SIZE) $(EXE)

# Library targets

lib:	$(OBJ) $(EXTRA_LINK_DEPENDS)

	$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)

shlib:	$(OBJ) $(EXTRA_LINK_DEPENDS)

	$(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o $(EXE) \

        $(OBJ) $(EXTRA_LIB) $(LIB)

# Compilation rules

%.o:%.cpp $(EXTRA_CPP_DEPENDS)

	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<

%.d:%.cpp $(EXTRA_CPP_DEPENDS)

	$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@

%.o:%.cu $(EXTRA_CPP_DEPENDS)

	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<

# Individual dependencies

depend : fastdep.exe $(SRC)

	@./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1

fastdep.exe: ../DEPEND/fastdep.c

	cc -O -o $@ $<

sinclude .depend

/* ----------------------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.

   ------------------------------------------------------------------------- */

/*  -----------------------------------------------------------------------

   Contributed by Stefan Paquay @ Brandeis University

   Thanks to Liesbeth Janssen @ Eindhoven University for useful discussions!

   ----------------------------------------------------------------------- */

#include <stdio.h>

#include <string.h>

#include "fix_active.h"

#include "atom.h"

#include "atom_vec_ellipsoid.h"

#include "citeme.h"

#include "comm.h"

#include "error.h"

#include "force.h"

#include "group.h"

#include "math.h"

#include "math_const.h"

#include "math_extra.h"

#include "math_vector.h"

#include "modify.h"

#include "random_mars.h"

#include "respa.h"

#include "update.h"

using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathConst;

/*

  Might be used later, who knows...

static const char* cite_fix_active =

  "fix active command:\n\n"

  "@article{paquay-2018,\n"

  "   author        = {Paquay, Stefan},\n"

  "   doi           = {},\n"

  "   issn          = {},\n"

  "   journal       = {},\n"

  "   month         = ,\n"

  "   number        = ,\n"

  "   pages         = ,\n"

  "   title         = ,\n"

  "   volume        = ,\n"

  "   year          = \n"

  "}\n\n";

*/

static constexpr const bool debug_out = false;

FixActive::FixActive( LAMMPS *lmp, int narg, char **argv )

  : Fix(lmp, narg, argv)

{

  // if( lmp->citeme) lmp->citeme->add(cite_fix_active);

  if( narg < 4 ) error->all(FLERR, "Illegal fix active command");

  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);

  if (!av) error->all(FLERR, "FixActive requires atom_style ellipsoid");

  if( debug_out && comm->me == 0 ){

    fprintf(screen, "This is fix active, narg is %d\n", narg );

    fprintf(screen, "args:");

    for( int i = 0; i < narg; ++i ){

      fprintf(screen, " %s", argv[i]);

    }

    fprintf(screen, "\n");

  }

  // args: fix ID all active magnitude prop1 prop2 prop3

  // Optional args are

  magnitude = force->numeric( FLERR, argv[3] );

}

FixActive::~FixActive()

{}

int FixActive::setmask()

{

  int mask = 0;

  mask |= POST_FORCE;

  return mask;

}

double FixActive::memory_usage()

{

  double bytes = 0.0;

  return bytes;

}

void FixActive::post_force(int /* vflag */ )

{

  // Then do the rest:

  double **f = atom->f;

  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  if (igroup == atom->firstgroup) nlocal = atom->nfirst;

  AtomVecEllipsoid::Bonus *bonus = av->bonus;

  // Add the active force to the atom force:

  for( int i = 0; i < nlocal; ++i ){

    if( mask[i] & groupbit ){

      double f_act[3] = { 0.0, 0.0, 1.0 };

      double f_rot[3];

      int* ellipsoid = atom->ellipsoid;

      double *quat  = bonus[ellipsoid[i]].quat;

      tagint *tag = atom->tag;

      double Q[3][3];

      MathExtra::quat_to_mat( quat, Q );

      MathExtra::matvec( Q, f_act, f_rot );

      if (debug_out && comm->me == 0) {

        // Magical reference particle:

        if (tag[i] == 12) {

          fprintf(screen, "rotation quaternion: (%f %f %f %f)\n",

                  quat[0], quat[1], quat[2], quat[3]);

          fprintf(screen, "rotation matrix: / %f %f %f \\\n",

                  Q[0][0] ,Q[0][1], Q[0][2]);

          fprintf(screen, "                 | %f %f %f |\n",

                  Q[1][0] ,Q[1][1], Q[1][2]);

          fprintf(screen, "                 \\ %f %f %f /\n",

                  Q[2][0] ,Q[2][1], Q[2][2]);

          fprintf(screen, "Active force on atom %d: (%f %f %f)\n",

                  tag[i], f_rot[0], f_rot[1], f_rot[2]);

          fprintf(screen, "  Total force before: (%f %f %f)\n",

                  f[i][0], f[i][1], f[i][2]);

        }

      }

      f[i][0] += magnitude * f_rot[0];

      f[i][1] += magnitude * f_rot[1];

      f[i][2] += magnitude * f_rot[2];

    }

  }

}