libghemical-dev 3.0.0-4.1build1 / usr / include / ghemical / moldyn.h

// MOLDYN.H : molecular dynamics classes.

// Copyright (C) 1998 Tommi Hassinen.

// This package is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.

// This package is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with this package; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

/*################################################################################################*/

#ifndef MOLDYN_H
#define MOLDYN_H

class moldyn_param;

class moldyn;

/*################################################################################################*/

#include "engine.h"
#include "eng1_sf.h"

#include <cstring>	// this is for std::strcpy ; 20071024

/*################################################################################################*/

class moldyn_param
{
	protected:
	
	bool confirm;		// for GUI only...
	
	public:
	
	i32s nsteps_h;		// heating
	i32s nsteps_e;		// equilibration
	i32s nsteps_s;		// simulation
	i32s nsteps_c;		// cooling
	
	f64 timestep;		// simulation timestep [fs]
	
	f64 target_T;		// target temperature [K]
	f64 T_rtime_hc;		// temperature relaxation time [fs] heating/cooling
	f64 T_rtime_es;		// temperature relaxation time [fs] equilibration/simulation
	
	f64 target_P;		// target pressure [bar]
	f64 P_rtime;		// pressure relaxation time [fs]
	f64 P_beta;		// isothermal compressibility [1/bar]
	
	bool constant_T;
	bool constant_P;
	
	char filename[256];
	
	public:
	
	moldyn_param(setup * su)
	{
		confirm = false;
		
	// set the defaults...
	// ^^^^^^^^^^^^^^^^^^^
		
		nsteps_h = 5000;
		nsteps_e = 5000;
		nsteps_s = 18000;
		nsteps_c = 2000;
		
		timestep = 0.5;
	setup1_sf * susf = dynamic_cast<setup1_sf *>(su);
	if (susf != NULL) timestep = 5.0;	// override...
		
		target_T = 300.0;
		T_rtime_hc = 10.0 * timestep;
		T_rtime_es = 100.0 * timestep;
		
		target_P = 1.000;
		P_rtime = 100.0 * timestep;
		P_beta = 4.57e-5;
		
		constant_T = true;
		constant_P = false;
		
		std::strcpy(filename, "untitled.traj");
	}
	
	~moldyn_param(void) { }
	
	bool GetConfirm(void) { return confirm; }
	void Confirm(void) { confirm = true; }
};

/*################################################################################################*/

/**	A molecular dynamics class.
	
	This is a "velocity-Verlet"-type integrator...
	Allen MP, Tildesley DJ : "##Computer Simulation of Liquids", Clarendon Press, 1987
	
	So far very long simulations have not been tested, so possile translation/rotation 
	problems are not solved... current solution is to start simulations from 0 K -> 
	no translation/rotation in the initial state -> conservation of linear and angular 
	momentum -> no need to worry at all... but for how long it will work???
*/

class moldyn
{
	protected:
	
	engine * eng;
	
	f64 * vel;			// [1.0e+3 m/s]
	f64 * acc;			// [1.0e+12 m/s^2]
	
	f64 * mass;			// [kg/mol]
	
	char * locked;
	int num_locked;
	
	f64 tstep1;		// timestep [fs]
	f64 tstep2;		// timestep ^ 2
	
	f64 ekin;
	f64 epot;
	
	i32s step_counter;
	
	f64 sum_of_masses;	// this is for density...
	
	friend void model::WriteTrajectoryFrame(ofstream &, moldyn *);
	
	public:
	
	f64 target_temperature;		// [K]
	f64 temperature_rtime;		// [fs]
	
	f64 target_pressure;		// [bar]
	f64 pressure_rtime;		// [fs]
	f64 isoth_compr;		// [1/bar]
	
	f64 saved_pressure;
	f64 saved_density;
	
	public:
	
	moldyn(engine *, f64);
	virtual ~moldyn(void);
	
	f64 GetEKin(void) { return ekin; }
	f64 GetEPot(void) { return epot; }
	
	virtual void TakeMDStep(bool, bool);
	
	f64 KineticEnergy(f64 * = NULL);
	
	f64 ConvTempEKin(f64);
	f64 ConvEKinTemp(f64);
	
	void SetEKin(f64);
};

/*################################################################################################*/

#endif	// MOLDYN_H

// eof