Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F75761293
fix_lambdah_calc.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Aug 4, 02:46
Size
32 KB
Mime Type
text/x-c
Expires
Tue, Aug 6, 02:46 (2 d)
Engine
blob
Format
Raw Data
Handle
19528654
Attached To
rLAMMPS lammps
fix_lambdah_calc.cpp
View Options
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Maziar Heidari (Max Planck Institute for Polymer Research)
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_lambdah_calc.h"
#include "atom.h"
#include "input.h"
#include "variable.h"
#include "domain.h"
#include "lattice.h"
#include "update.h"
#include "modify.h"
#include "output.h"
#include "respa.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "comm.h"
#include "citeme.h"
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
using
namespace
MathConst
;
#define BIG MAXTAGINT
static
const
char
cite_HAdResS
[]
=
"@Article{Heidari et al.2016
\n
"
" author = {M. Heidari, R. Cortes-Huerto, D. Donadio and R. Potestio},
\n
"
" title = {Accurate and general treatment of electrostatic interaction in Hamiltonian adaptive resolution simulations},
\n
"
" journal = {Eur. Phys. J. Special Topics},
\n
"
" year = 2016,
\n
"
" volume = Submitted,
\n
"
" pages = {}
\n
"
"}
\n\n
"
;
/* ---------------------------------------------------------------------- */
FixLambdaHCalc
::
FixLambdaHCalc
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
)
{
if
(
lmp
->
citeme
)
lmp
->
citeme
->
add
(
cite_HAdResS
);
massproc_H
=
NULL
;
masstotal_H
=
NULL
;
com_H
=
NULL
;
comall_H
=
NULL
;
molmap_H
=
NULL
;
lambdaCM
=
NULL
;
gradlambdaCM
=
NULL
;
Comp_Density_Num_H
=
NULL
;
Comp_Density_Num_all_H
=
NULL
;
Mean_Density_H
=
NULL
;
Int_Mean_Density_H
=
NULL
;
Mean_Comp_Density_Conv_H
=
NULL
;
grad_Comp_Density_Conv_H
=
NULL
;
Mean_grad_Comp_Density_Conv_H
=
NULL
;
me
=
comm
->
me
;
if
(
narg
<
8
)
error
->
all
(
FLERR
,
"Illegal fix lambdah/calc command"
);
atom
->
nmoltypesH
=
force
->
numeric
(
FLERR
,
arg
[
3
]);
Length_Hyb
=
force
->
numeric
(
FLERR
,
arg
[
4
]);
Length_ATRegion
=
force
->
numeric
(
FLERR
,
arg
[
5
]);
Pressure_Comp_Flag
=
force
->
numeric
(
FLERR
,
arg
[
6
]);
Pressure_lambda_Increment
=
force
->
numeric
(
FLERR
,
arg
[
7
]);
Pressure_Update_Frequency
=
force
->
numeric
(
FLERR
,
arg
[
8
]);
Pressure_Update_Time_Begin
=
force
->
numeric
(
FLERR
,
arg
[
9
]);
Pressure_Update_Time_End
=
force
->
numeric
(
FLERR
,
arg
[
10
]);
if
(
strcmp
(
arg
[
11
],
"slab"
)
==
0
)
Hybrid_Style
=
0
;
else
if
(
strcmp
(
arg
[
11
],
"sphere"
)
==
0
)
Hybrid_Style
=
1
;
else
if
(
strcmp
(
arg
[
11
],
"cylinder"
)
==
0
)
Hybrid_Style
=
2
;
else
error
->
all
(
FLERR
,
"Illegal fix lambdah/calc command"
);
Density_Comp_Flag
=
force
->
numeric
(
FLERR
,
arg
[
12
]);
Density_Bin_Size
=
force
->
numeric
(
FLERR
,
arg
[
13
]);
Density_Update_Frequency
=
force
->
numeric
(
FLERR
,
arg
[
14
]);
Density_Update_Time_Begin
=
force
->
numeric
(
FLERR
,
arg
[
15
]);
Density_Update_Time_End
=
force
->
numeric
(
FLERR
,
arg
[
16
]);
Density_Sigma_Gauss
=
force
->
numeric
(
FLERR
,
arg
[
17
]);
Density_Gauss_Int_Range
=
force
->
numeric
(
FLERR
,
arg
[
18
]);
Density_Ref
=
force
->
numeric
(
FLERR
,
arg
[
19
]);
Comp_Density_Scaling_factor_H
=
force
->
numeric
(
FLERR
,
arg
[
20
]);
Load_File_Flag
=
force
->
numeric
(
FLERR
,
arg
[
21
]);
x0lo
=
domain
->
boxlo
[
0
];
x0hi
=
domain
->
boxhi
[
0
];
x1lo
=
domain
->
boxlo
[
1
];
x1hi
=
domain
->
boxhi
[
1
];
x2lo
=
domain
->
boxlo
[
2
];
x2hi
=
domain
->
boxhi
[
2
];
center_box
[
0
]
=
(
x0hi
+
x0lo
)
/
2.0
;
center_box
[
1
]
=
(
x1hi
+
x1lo
)
/
2.0
;
center_box
[
2
]
=
(
x2hi
+
x2lo
)
/
2.0
;
x0BoxSize
=
x0hi
-
x0lo
;
Length_CGRegion
=
x0BoxSize
-
2
*
Length_Hyb
-
Length_ATRegion
;
R_Start_Hybrid_1
=
x0lo
+
Length_CGRegion
/
2.0
;
R_Start_Hybrid_2
=
x0lo
+
x0BoxSize
-
(
Length_CGRegion
/
2.0
+
Length_Hyb
);
S_Start_Hybrid
=
Length_ATRegion
;
Pressure_Bin_Num
=
1.0
/
Pressure_lambda_Increment
;
xmin_AT
=
R_Start_Hybrid_1
+
Length_Hyb
;
xmax_AT
=
R_Start_Hybrid_1
+
Length_Hyb
+
Length_ATRegion
;
if
(
Hybrid_Style
==
0
)
Density_Bin_Num
=
floor
(
x0BoxSize
/
Density_Bin_Size
);
else
if
(
Hybrid_Style
==
1
)
Density_Bin_Num
=
floor
(
sqrt
(
pow
(
0.5
*
(
x0hi
-
x0lo
),
2.0
)
+
pow
(
0.5
*
(
x1hi
-
x1lo
),
2.0
)
+
pow
(
0.5
*
(
x2hi
-
x2lo
),
2.0
))
/
Density_Bin_Size
);
else
if
(
Hybrid_Style
==
2
)
Density_Bin_Num
=
floor
(
sqrt
(
pow
(
0.5
*
(
x0hi
-
x0lo
),
2.0
)
+
pow
(
0.5
*
(
x1hi
-
x1lo
),
2.0
))
/
Density_Bin_Size
);
Comp_Counter_H
=
0
;
Density_Counter_H
=
0
;
Density_Compensation_Run
=
0
;
memory
->
create
(
Comp_Density_Num_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:Comp_Density_Num_H"
);
memory
->
create
(
Comp_Density_Num_all_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:Comp_Density_Num_all_H"
);
if
(
me
==
0
){
if
(
screen
){
fprintf
(
screen
,
"nmoltypes= %d
\n
"
,
atom
->
nmoltypesH
);
fprintf
(
screen
,
"Length_Hyb= %f
\n
"
,
Length_Hyb
);
fprintf
(
screen
,
"Length_ATRegion= %f
\n
"
,
Length_ATRegion
);
fprintf
(
screen
,
"Pressure_Comp_Flag= %d
\n
"
,
Pressure_Comp_Flag
);
fprintf
(
screen
,
"Pressure_lambda_Increment= %f
\n
"
,
Pressure_lambda_Increment
);
fprintf
(
screen
,
"Pressure_Update_Frequency= %d
\n
"
,
Pressure_Update_Frequency
);
fprintf
(
screen
,
"Pressure_Update_Time_Begin= %d
\n
"
,
Pressure_Update_Time_Begin
);
fprintf
(
screen
,
"Pressure_Update_Time_End= %d
\n
"
,
Pressure_Update_Time_End
);
fprintf
(
screen
,
"Density_Comp_Flag= %d
\n
"
,
Density_Comp_Flag
);
fprintf
(
screen
,
"Density_Bin_Size= %f
\n
"
,
Density_Bin_Size
);
fprintf
(
screen
,
"Density_Update_Frequency= %d
\n
"
,
Density_Update_Frequency
);
fprintf
(
screen
,
"Density_Update_Time_Begin= %d
\n
"
,
Density_Update_Time_Begin
);
fprintf
(
screen
,
"Density_Update_Time_End= %d
\n
"
,
Density_Update_Time_End
);
fprintf
(
screen
,
"Density_Sigma_Gauss= %f
\n
"
,
Density_Sigma_Gauss
);
fprintf
(
screen
,
"Density_Gauss_Int_Range= %d
\n
"
,
Density_Gauss_Int_Range
);
fprintf
(
screen
,
"Density_Ref= %f
\n
"
,
Density_Ref
);
fprintf
(
screen
,
"Comp_Density_Scaling_factor_H = %f
\n
"
,
Comp_Density_Scaling_factor_H
);
fprintf
(
screen
,
"Load_File_Flag = %d
\n
"
,
Load_File_Flag
);
fprintf
(
screen
,
"Center_box = %f %f %f
\n
"
,
center_box
[
0
],
center_box
[
1
],
center_box
[
2
]);
fprintf
(
screen
,
"Density_Bin_Size = %f
\n
"
,
Density_Bin_Size
);
fprintf
(
screen
,
"Density_Bin_Num = %d
\n
"
,
Density_Bin_Num
);
fprintf
(
screen
,
"x0lo= %f
\n
"
,
x0lo
);
fprintf
(
screen
,
"x0hi= %f
\n
"
,
x0hi
);
fprintf
(
screen
,
"x0BoxSize= %f
\n
"
,
x0BoxSize
);
fprintf
(
screen
,
"d1= %f
\n
"
,
R_Start_Hybrid_1
);
fprintf
(
screen
,
"d2= %f
\n
"
,
R_Start_Hybrid_2
);
fprintf
(
screen
,
"moltype%d
\n
"
,
atom
->
nmoltypesH
);
}
if
(
logfile
){
fprintf
(
logfile
,
"nmoltypes= %d
\n
"
,
atom
->
nmoltypesH
);
fprintf
(
logfile
,
"Length_Hyb= %f
\n
"
,
Length_Hyb
);
fprintf
(
logfile
,
"Length_ATRegion= %f
\n
"
,
Length_ATRegion
);
fprintf
(
logfile
,
"Pressure_Comp_Flag= %d
\n
"
,
Pressure_Comp_Flag
);
fprintf
(
logfile
,
"Pressure_lambda_Increment= %f
\n
"
,
Pressure_lambda_Increment
);
fprintf
(
logfile
,
"Pressure_Update_Frequency= %d
\n
"
,
Pressure_Update_Frequency
);
fprintf
(
logfile
,
"Pressure_Update_Time_Begin= %d
\n
"
,
Pressure_Update_Time_Begin
);
fprintf
(
logfile
,
"Pressure_Update_Time_End= %d
\n
"
,
Pressure_Update_Time_End
);
fprintf
(
logfile
,
"Density_Comp_Flag= %d
\n
"
,
Density_Comp_Flag
);
fprintf
(
logfile
,
"Density_Bin_Size= %f
\n
"
,
Density_Bin_Size
);
fprintf
(
logfile
,
"Density_Update_Frequency= %d
\n
"
,
Density_Update_Frequency
);
fprintf
(
logfile
,
"Density_Update_Time_Begin= %d
\n
"
,
Density_Update_Time_Begin
);
fprintf
(
logfile
,
"Density_Update_Time_End= %d
\n
"
,
Density_Update_Time_End
);
fprintf
(
logfile
,
"Density_Sigma_Gauss= %f
\n
"
,
Density_Sigma_Gauss
);
fprintf
(
logfile
,
"Density_Gauss_Int_Range= %d
\n
"
,
Density_Gauss_Int_Range
);
fprintf
(
logfile
,
"Density_Ref= %f
\n
"
,
Density_Ref
);
fprintf
(
logfile
,
"Comp_Density_Scaling_factor_H = %f
\n
"
,
Comp_Density_Scaling_factor_H
);
fprintf
(
logfile
,
"Load_File_Flag = %d
\n
"
,
Load_File_Flag
);
fprintf
(
logfile
,
"Center_box = %f %f %f
\n
"
,
center_box
[
0
],
center_box
[
1
],
center_box
[
2
]);
fprintf
(
logfile
,
"Density_Bin_Size = %f
\n
"
,
Density_Bin_Size
);
fprintf
(
logfile
,
"Density_Bin_Num = %d
\n
"
,
Density_Bin_Num
);
fprintf
(
logfile
,
"x0lo= %f
\n
"
,
x0lo
);
fprintf
(
logfile
,
"x0hi= %f
\n
"
,
x0hi
);
fprintf
(
logfile
,
"x0BoxSize= %f
\n
"
,
x0BoxSize
);
fprintf
(
logfile
,
"d1= %f
\n
"
,
R_Start_Hybrid_1
);
fprintf
(
logfile
,
"d2= %f
\n
"
,
R_Start_Hybrid_2
);
fprintf
(
logfile
,
"moltype%d
\n
"
,
atom
->
nmoltypesH
);
}
}
memory
->
create
(
Int_Mean_Density_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:Int_Mean_Density_H"
);
memory
->
create
(
Mean_Density_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:Mean_Density_H"
);
memory
->
create
(
Mean_Comp_Density_Conv_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:Mean_Comp_Density_Conv_H"
);
memory
->
create
(
grad_Comp_Density_Conv_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:grad_Comp_Density_Conv_H"
);
memory
->
create
(
Mean_grad_Comp_Density_Conv_H
,
Density_Bin_Num
,
atom
->
nmoltypesH
+
1
,
"lambdaH/calc:Mean_grad_Comp_Density_Conv_H"
);
int
This_Step
=
update
->
ntimestep
;
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
){
for
(
int
j
=
0
;
j
<
atom
->
nmoltypesH
;
j
++
){
Comp_Density_Num_H
[
i
][
j
]
=
0
;
Comp_Density_Num_all_H
[
i
][
j
]
=
0
;
Int_Mean_Density_H
[
i
][
j
]
=
0
;
Mean_Density_H
[
i
][
j
]
=
0
;
Mean_Comp_Density_Conv_H
[
i
][
j
]
=
0
;
Mean_grad_Comp_Density_Conv_H
[
i
][
j
]
=
0
;
grad_Comp_Density_Conv_H
[
i
][
j
]
=
0
;
}
}
if
((
This_Step
>=
Density_Update_Time_End
||
Load_File_Flag
)
&&
Density_Comp_Flag
!=
0
)
Load_Compensation_Density
();
if
(
atom
->
molecular
==
0
)
error
->
all
(
FLERR
,
"Compute com/molecule requires molecular atom style"
);
array_flag
=
1
;
size_array_cols
=
3
;
extarray
=
0
;
// setup molecule-based data
nmolecules
=
molecules_in_group
(
idlo
,
idhi
);
size_array_rows
=
nmolecules
;
//printf("Nmolecules= %d\n",nmolecules);
memory
->
create
(
massproc_H
,
nmolecules
,
"lambdaH/calc:massproc_H"
);
memory
->
create
(
masstotal_H
,
nmolecules
,
"lambdaH/calc:masstotal_H"
);
memory
->
create
(
com_H
,
nmolecules
,
3
,
"lambdaH/calc:com_H"
);
memory
->
create
(
comall_H
,
nmolecules
,
3
,
"lambdaH/calc:comall_H"
);
memory
->
create
(
lambdaCM
,
nmolecules
,
"lambdaH/calc:lambdaCM"
);
memory
->
create
(
gradlambdaCM
,
nmolecules
,
3
,
"lambdaH/calc:gradlambdaCM"
);
// compute masstotal for each molecule
int
*
mask
=
atom
->
mask
;
tagint
*
molecule
=
atom
->
molecule
;
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
tagint
imol
;
double
massone
;
for
(
int
i
=
0
;
i
<
nmolecules
;
i
++
)
massproc_H
[
i
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
{
if
(
rmass
)
massone
=
rmass
[
i
];
else
massone
=
mass
[
type
[
i
]];
imol
=
molecule
[
i
];
if
(
molmap_H
)
imol
=
molmap_H
[
imol
-
idlo
];
else
imol
--
;
massproc_H
[
imol
]
+=
massone
;
}
}
MPI_Allreduce
(
massproc_H
,
masstotal_H
,
nmolecules
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
int
*
replambdaH
=
atom
->
replambdaH
;
int
pass
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
replambdaH
[
i
]
!=
0
)
pass
=
1
;
if
(
pass
==
0
)
error
->
all
(
FLERR
,
"No representative atom (replambdaH) has been defined!"
);
}
/* ---------------------------------------------------------------------- */
FixLambdaHCalc
::~
FixLambdaHCalc
()
{
memory
->
destroy
(
massproc_H
);
memory
->
destroy
(
masstotal_H
);
memory
->
destroy
(
com_H
);
memory
->
destroy
(
comall_H
);
memory
->
destroy
(
molmap_H
);
memory
->
destroy
(
lambdaCM
);
memory
->
destroy
(
gradlambdaCM
);
memory
->
destroy
(
Comp_Density_Num_H
);
memory
->
destroy
(
Comp_Density_Num_all_H
);
memory
->
destroy
(
Mean_Density_H
);
memory
->
destroy
(
Int_Mean_Density_H
);
memory
->
destroy
(
Mean_Comp_Density_Conv_H
);
memory
->
destroy
(
grad_Comp_Density_Conv_H
);
memory
->
destroy
(
Mean_grad_Comp_Density_Conv_H
);
}
/* ---------------------------------------------------------------------- */
int
FixLambdaHCalc
::
setmask
()
{
int
mask
=
0
;
//mask |= PRE_FORCE;
// mask |= PRE_NEIGHBOR;
mask
|=
POST_INTEGRATE
;
mask
|=
THERMO_ENERGY
;
// mask |= POST_FORCE_RESPA;
mask
|=
MIN_PRE_FORCE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixLambdaHCalc
::
init
()
{
int
ntmp
=
molecules_in_group
(
idlo
,
idhi
);
if
(
ntmp
!=
nmolecules
)
error
->
all
(
FLERR
,
"Molecule count changed in compute com/molecule"
);
}
/* ---------------------------------------------------------------------- */
void
FixLambdaHCalc
::
setup
(
int
vflag
)
{
// if (strstr(update->integrate_style,"verlet"))
//pre_force(vflag);
post_integrate
();
// else {
// ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
// post_force_respa(vflag,nlevels_respa-1,0);
// ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
// }
}
/* ---------------------------------------------------------------------- */
void
FixLambdaHCalc
::
post_integrate
()
{
tagint
imol
;
double
massone
;
double
unwrap
[
3
];
//invoked_array = update->ntimestep;
for
(
int
i
=
0
;
i
<
nmolecules
;
i
++
)
com_H
[
i
][
0
]
=
com_H
[
i
][
1
]
=
com_H
[
i
][
2
]
=
0.0
;
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
tagint
*
molecule
=
atom
->
molecule
;
int
*
type
=
atom
->
type
;
imageint
*
image
=
atom
->
image
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
*
lambdaH
=
atom
->
lambdaH
;
double
**
gradlambdaH
=
atom
->
gradlambdaH
;
double
**
comH
=
atom
->
comH
;
int
This_Step
=
update
->
ntimestep
;
if
(
This_Step
>=
Density_Update_Time_Begin
&&
This_Step
<
Density_Update_Time_End
&&
Density_Comp_Flag
==
1
){
Density_Compensation_Run
=
1
;
if
(
me
==
0
&&
This_Step
==
Density_Update_Time_Begin
){
if
(
screen
)
fprintf
(
screen
,
"
\n
Start of constant-density route
\n
"
);
if
(
logfile
)
fprintf
(
logfile
,
"
\n
Start of constant-density route
\n
"
);
Clear_File_Compensation_Density
();
}
}
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
if
(
rmass
)
massone
=
rmass
[
i
];
else
massone
=
mass
[
type
[
i
]];
imol
=
molecule
[
i
];
if
(
molmap_H
)
imol
=
molmap_H
[
imol
-
idlo
];
else
imol
--
;
domain
->
unmap
(
x
[
i
],
image
[
i
],
unwrap
);
massone
/=
masstotal_H
[
imol
];
com_H
[
imol
][
0
]
+=
unwrap
[
0
]
*
massone
;
com_H
[
imol
][
1
]
+=
unwrap
[
1
]
*
massone
;
com_H
[
imol
][
2
]
+=
unwrap
[
2
]
*
massone
;
}
MPI_Allreduce
(
&
com_H
[
0
][
0
],
&
comall_H
[
0
][
0
],
3
*
nmolecules
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
double
xtmp
,
ytmp
,
ztmp
;
double
delx
,
dely
,
delz
;
double
Grad_Factor
;
double
rdiff
;
for
(
int
i
=
0
;
i
<
nmolecules
;
i
++
){
domain
->
remap
(
comall_H
[
i
]);
if
(
Hybrid_Style
==
0
){
xtmp
=
comall_H
[
i
][
0
];
if
(
xtmp
<
R_Start_Hybrid_1
){
lambdaCM
[
i
]
=
0.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
xtmp
>=
R_Start_Hybrid_1
&&
xtmp
<
R_Start_Hybrid_1
+
Length_Hyb
){
cosx
=
cos
(
MY_PI
*
(
xtmp
-
R_Start_Hybrid_1
+
Length_Hyb
)
/
(
2.0
*
Length_Hyb
));
sinx
=
sin
(
MY_PI
*
(
xtmp
-
R_Start_Hybrid_1
+
Length_Hyb
)
/
(
2.0
*
Length_Hyb
));
lambdaCM
[
i
]
=
cosx
*
cosx
;
gradlambdaCM
[
i
][
0
]
=
-
(
MY_PI
/
Length_Hyb
)
*
sinx
*
cosx
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
xtmp
>=
R_Start_Hybrid_1
+
Length_Hyb
&&
xtmp
<
R_Start_Hybrid_2
){
lambdaCM
[
i
]
=
1.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
xtmp
>=
R_Start_Hybrid_2
&&
xtmp
<
R_Start_Hybrid_2
+
Length_Hyb
){
cosx
=
cos
(
MY_PI
*
(
xtmp
-
R_Start_Hybrid_2
)
/
(
2.0
*
Length_Hyb
));
sinx
=
sin
(
MY_PI
*
(
xtmp
-
R_Start_Hybrid_2
)
/
(
2.0
*
Length_Hyb
));
lambdaCM
[
i
]
=
cosx
*
cosx
;
gradlambdaCM
[
i
][
0
]
=
-
(
MY_PI
/
Length_Hyb
)
*
sinx
*
cosx
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
xtmp
>=
R_Start_Hybrid_2
+
Length_Hyb
){
lambdaCM
[
i
]
=
0.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
}
else
if
(
Hybrid_Style
==
1
){
xtmp
=
comall_H
[
i
][
0
];
ytmp
=
comall_H
[
i
][
1
];
ztmp
=
comall_H
[
i
][
2
];
delx
=
(
xtmp
-
center_box
[
0
]);
dely
=
(
ytmp
-
center_box
[
1
]);
delz
=
(
ztmp
-
center_box
[
2
]);
r
=
sqrt
(
delx
*
delx
+
dely
*
dely
+
delz
*
delz
);
if
(
r
<
S_Start_Hybrid
){
lambdaCM
[
i
]
=
1.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
r
>=
S_Start_Hybrid
&&
r
<
S_Start_Hybrid
+
Length_Hyb
)
{
rdiff
=
MY_PI
*
(
r
-
S_Start_Hybrid
)
/
(
2.0
*
Length_Hyb
);
cosr
=
cos
(
rdiff
);
sinr
=
sin
(
rdiff
);
lambdaCM
[
i
]
=
cosr
*
cosr
;
Grad_Factor
=
-
MY_PI
*
sinr
*
cosr
/
(
Length_Hyb
*
r
);
gradlambdaCM
[
i
][
0
]
=
Grad_Factor
*
delx
;
gradlambdaCM
[
i
][
1
]
=
Grad_Factor
*
dely
;
gradlambdaCM
[
i
][
2
]
=
Grad_Factor
*
delz
;
}
else
if
(
r
>=
S_Start_Hybrid
+
Length_Hyb
){
lambdaCM
[
i
]
=
0.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
}
else
if
(
Hybrid_Style
==
2
){
xtmp
=
comall_H
[
i
][
0
];
ytmp
=
comall_H
[
i
][
1
];
delx
=
(
xtmp
-
center_box
[
0
]);
dely
=
(
ytmp
-
center_box
[
1
]);
r
=
sqrt
(
delx
*
delx
+
dely
*
dely
);
if
(
r
<
S_Start_Hybrid
){
lambdaCM
[
i
]
=
1.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
r
>=
S_Start_Hybrid
&&
r
<
S_Start_Hybrid
+
Length_Hyb
)
{
rdiff
=
MY_PI
*
(
r
-
S_Start_Hybrid
)
/
(
2.0
*
Length_Hyb
);
cosr
=
cos
(
rdiff
);
sinr
=
sin
(
rdiff
);
lambdaCM
[
i
]
=
cosr
*
cosr
;
Grad_Factor
=
-
MY_PI
*
sinr
*
cosr
/
(
Length_Hyb
*
r
);
gradlambdaCM
[
i
][
0
]
=
Grad_Factor
*
delx
;
gradlambdaCM
[
i
][
1
]
=
Grad_Factor
*
dely
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
else
if
(
r
>=
S_Start_Hybrid
+
Length_Hyb
){
lambdaCM
[
i
]
=
0.0
;
gradlambdaCM
[
i
][
0
]
=
0.0
;
gradlambdaCM
[
i
][
1
]
=
0.0
;
gradlambdaCM
[
i
][
2
]
=
0.0
;
}
}
}
int
ibin
,
imoltypeH
;
int
*
moltypeH
=
atom
->
moltypeH
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
imol
=
molecule
[
i
];
imoltypeH
=
moltypeH
[
i
];
if
(
molmap_H
)
imol
=
molmap_H
[
imol
-
idlo
];
else
imol
--
;
lambdaH
[
i
]
=
lambdaCM
[
imol
];
gradlambdaH
[
i
][
0
]
=
gradlambdaCM
[
imol
][
0
];
gradlambdaH
[
i
][
1
]
=
gradlambdaCM
[
imol
][
1
];
gradlambdaH
[
i
][
2
]
=
gradlambdaCM
[
imol
][
2
];
// if(gradlambdaH[i][2]!=0) printf("gradlambdaH[i][2] = %f\n",gradlambdaH[i][2]);
// gradlambdaH[i][2] = 0.0;
//gradlambdaH[i] = 0;
// if(replambdaH[i] == 1){
comH
[
i
][
0
]
=
comall_H
[
imol
][
0
];
comH
[
i
][
1
]
=
comall_H
[
imol
][
1
];
comH
[
i
][
2
]
=
comall_H
[
imol
][
2
];
// }
if
(
Density_Comp_Flag
!=
0
&&
Density_Compensation_Run
!=
0
){
if
(
Hybrid_Style
==
0
){
// xtmp = comall_H[imol][0]-x0lo;
xtmp
=
comH
[
i
][
0
]
-
x0lo
;
// if(xtmp <= 0.5*x0BoxSize)ibin = floor(xtmp/Density_Bin_Size);
// else ibin = floor((x0BoxSize-xtmp)/Density_Bin_Size);
ibin
=
floor
(
xtmp
/
Density_Bin_Size
);
}
else
if
(
Hybrid_Style
==
1
){
xtmp
=
comall_H
[
imol
][
0
];
ytmp
=
comall_H
[
imol
][
1
];
ztmp
=
comall_H
[
imol
][
2
];
delx
=
(
xtmp
-
center_box
[
0
]);
dely
=
(
ytmp
-
center_box
[
1
]);
delz
=
(
ztmp
-
center_box
[
2
]);
r
=
sqrt
(
delx
*
delx
+
dely
*
dely
+
delz
*
delz
);
ibin
=
floor
(
r
/
Density_Bin_Size
);
}
else
if
(
Hybrid_Style
==
2
){
xtmp
=
comall_H
[
imol
][
0
];
ytmp
=
comall_H
[
imol
][
1
];
delx
=
(
xtmp
-
center_box
[
0
]);
dely
=
(
ytmp
-
center_box
[
1
]);
r
=
sqrt
(
delx
*
delx
+
dely
*
dely
);
ibin
=
floor
(
r
/
Density_Bin_Size
);
}
if
(
ibin
>=
Density_Bin_Num
)
ibin
=
Density_Bin_Num
-
1
;
Comp_Density_Num_H
[
ibin
][
imoltypeH
-
1
]
++
;
}
}
// int Middle_Bin_Num = (Density_Bin_Num-1) / 2;
// if(Hybrid_Style == 0)for(int i = Middle_Bin_Num+1;i < Density_Bin_Num; i++)Comp_Density_Num_H[i] = Comp_Density_Num_H[Density_Bin_Num-1-i];
if
(
Density_Compensation_Run
)
Density_Counter_H
++
;
double
Density_Bin_Vol
,
exponent
;
double
Normalization
;
int
jmin
,
jmax
,
jj
;
if
(
This_Step
%
Density_Update_Frequency
==
0
&&
This_Step
>
Density_Update_Time_Begin
&&
Density_Comp_Flag
!=
0
&&
Density_Compensation_Run
!=
0
){
// if(atom->nmoltypesH ==1) MPI_Allreduce(&Comp_Density_Num_H[0][0],&Comp_Density_Num_all_H[0][0],Density_Bin_Num,MPI_INT,MPI_SUM,world);
MPI_Allreduce
(
&
Comp_Density_Num_H
[
0
][
0
],
&
Comp_Density_Num_all_H
[
0
][
0
],
Density_Bin_Num
*
(
atom
->
nmoltypesH
+
1
),
MPI_INT
,
MPI_SUM
,
world
);
if
(
Hybrid_Style
==
0
)
Density_Bin_Vol
=
1.0
*
Density_Bin_Size
*
(
x1hi
-
x1lo
)
*
(
x2hi
-
x2lo
);
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
){
if
(
Hybrid_Style
==
1
)
Density_Bin_Vol
=
(
4.0
/
3.0
)
*
MY_PI
*
(
pow
(((
i
+
1
)
*
Density_Bin_Size
),
3.0
)
-
pow
((
i
*
Density_Bin_Size
),
3.0
));
if
(
Hybrid_Style
==
2
)
Density_Bin_Vol
=
4.0
*
MY_PI
*
(
x2hi
-
x2lo
)
*
(
pow
(((
i
+
1
)
*
Density_Bin_Size
),
2.0
)
-
pow
((
i
*
Density_Bin_Size
),
2.0
));
for
(
int
j
=
0
;
j
<
atom
->
nmoltypesH
;
j
++
){
Mean_Density_H
[
i
][
j
]
=
Comp_Density_Num_all_H
[
i
][
j
]
/
(
Density_Counter_H
*
Density_Bin_Vol
);
}
}
Density_Counter_H
=
0
;
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
){
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
){
Normalization
=
0
;
jmin
=
i
-
Density_Gauss_Int_Range
*
floor
(
Density_Sigma_Gauss
/
Density_Bin_Size
);
jmax
=
i
+
Density_Gauss_Int_Range
*
floor
(
Density_Sigma_Gauss
/
Density_Bin_Size
);
Mean_Comp_Density_Conv_H
[
i
][
k
]
=
0
;
if
(
Hybrid_Style
!=
0
){
if
(
jmin
<
0
)
jmin
=
0
;
if
(
jmax
>=
Density_Bin_Num
)
jmax
=
Density_Bin_Num
-
1
;
}
for
(
int
j
=
jmin
;
j
<=
jmax
;
j
++
){
jj
=
j
;
if
(
Hybrid_Style
==
0
){
if
(
j
<
0
)
jj
=
Density_Bin_Num
+
j
;
if
(
j
>=
Density_Bin_Num
)
jj
=
j
-
Density_Bin_Num
;
}
exponent
=
(
i
-
j
)
*
Density_Bin_Size
/
Density_Sigma_Gauss
;
exponent
=
pow
(
exponent
,
2.0
);
Normalization
+=
exp
(
-
exponent
)
*
Density_Bin_Size
;
Mean_Comp_Density_Conv_H
[
i
][
k
]
+=
Mean_Density_H
[
jj
][
k
]
*
exp
(
-
exponent
)
*
Density_Bin_Size
;
// Mean_Comp_Density_Conv_H[i] += Mean_Density_H[j] * exp(-exponent);
}
Mean_Comp_Density_Conv_H
[
i
][
k
]
/=
Normalization
;
}
}
Density_Fluctuation
=
0
;
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
)
Density_Fluctuation
+=
(
Density_Bin_Size
/
x0BoxSize
)
*
pow
((
Mean_Density_H
[
i
][
k
]
-
Density_Ref
)
/
Density_Ref
,
2.0
);
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
for
(
int
i
=
1
;
i
<
Density_Bin_Num
-
1
;
i
++
)
grad_Comp_Density_Conv_H
[
i
][
k
]
=
Comp_Density_Scaling_factor_H
*
(
Mean_Comp_Density_Conv_H
[
i
+
1
][
k
]
-
Mean_Comp_Density_Conv_H
[
i
-
1
][
k
])
/
(
2
*
Density_Bin_Size
*
Density_Ref
);
//for(int i = 1;i < Density_Bin_Num-1; i++)Mean_grad_Comp_Density_Conv_H[i] = (Comp_Counter_H * Mean_grad_Comp_Density_Conv_H[i] + grad_Comp_Density_Conv_H[i]) / (Comp_Counter_H + 1);
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
for
(
int
i
=
1
;
i
<
Density_Bin_Num
-
1
;
i
++
)
Mean_grad_Comp_Density_Conv_H
[
i
][
k
]
+=
grad_Comp_Density_Conv_H
[
i
][
k
];
// for(int i = Middle_Bin_Num+1;i < Density_Bin_Num; i++)Comp_Density_Num_H[i] = Comp_Density_Num_H[Density_Bin_Num-1-i];
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
)
Int_Mean_Density_H
[
i
][
k
]
=
0
;
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
)
for
(
int
j
=
0
;
j
<=
i
;
j
++
)
Int_Mean_Density_H
[
i
][
k
]
+=
Mean_grad_Comp_Density_Conv_H
[
j
][
k
]
*
Density_Bin_Size
;
Comp_Counter_H
+=
1
;
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
){
Comp_Density_Num_H
[
i
][
k
]
=
0
;
Comp_Density_Num_all_H
[
i
][
k
]
=
0
;
}
}
if
(
This_Step
==
Density_Update_Time_End
&&
Density_Compensation_Run
!=
0
){
Density_Compensation_Run
=
0
;
if
(
me
==
0
){
if
(
screen
)
fprintf
(
screen
,
"
\n
End of constant-density route
\n
"
);
if
(
logfile
)
fprintf
(
logfile
,
"
\n
End of constant-density route
\n
"
);
}
}
if
(
me
==
0
&&
This_Step
%
Density_Update_Frequency
==
0
&&
Density_Compensation_Run
!=
0
)
Print_Compensation_Density
();
}
double
FixLambdaHCalc
::
memory_usage
()
{
double
bytes
=
(
bigint
)
nmolecules
*
2
*
sizeof
(
double
);
if
(
molmap_H
)
bytes
+=
(
idhi
-
idlo
+
1
)
*
sizeof
(
int
);
bytes
+=
(
bigint
)
nmolecules
*
2
*
3
*
sizeof
(
double
);
bytes
+=
(
bigint
)
nmolecules
*
2
*
3
*
sizeof
(
double
);
return
bytes
;
}
/* ---------------------------------------------------------------------- */
int
FixLambdaHCalc
::
molecules_in_group
(
tagint
&
idlo
,
tagint
&
idhi
)
{
int
i
;
memory
->
destroy
(
molmap_H
);
molmap_H
=
NULL
;
// find lo/hi molecule ID for any atom in group
// warn if atom in group has ID = 0
tagint
*
molecule
=
atom
->
molecule
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
tagint
lo
=
BIG
;
tagint
hi
=
-
BIG
;
int
flag
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
{
if
(
molecule
[
i
]
==
0
)
flag
=
1
;
lo
=
MIN
(
lo
,
molecule
[
i
]);
hi
=
MAX
(
hi
,
molecule
[
i
]);
}
}
int
flagall
;
MPI_Allreduce
(
&
flag
,
&
flagall
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flagall
&&
comm
->
me
==
0
)
error
->
warning
(
FLERR
,
"Atom with molecule ID = 0 included in "
"compute molecule group"
);
MPI_Allreduce
(
&
lo
,
&
idlo
,
1
,
MPI_LMP_TAGINT
,
MPI_MIN
,
world
);
MPI_Allreduce
(
&
hi
,
&
idhi
,
1
,
MPI_LMP_TAGINT
,
MPI_MAX
,
world
);
if
(
idlo
==
BIG
)
return
0
;
// molmap = vector of length nlen
// set to 1 for IDs that appear in group across all procs, else 0
tagint
nlen_tag
=
idhi
-
idlo
+
1
;
if
(
nlen_tag
>
MAXSMALLINT
)
error
->
all
(
FLERR
,
"Too many molecules for compute"
);
int
nlen
=
(
int
)
nlen_tag
;
memory
->
create
(
molmap_H
,
nlen
,
"lambdaH/calc:molmap"
);
for
(
i
=
0
;
i
<
nlen
;
i
++
)
molmap_H
[
i
]
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
molmap_H
[
molecule
[
i
]
-
idlo
]
=
1
;
int
*
molmapall
;
memory
->
create
(
molmapall
,
nlen
,
"compute:molmapall"
);
MPI_Allreduce
(
molmap_H
,
molmapall
,
nlen
,
MPI_INT
,
MPI_MAX
,
world
);
// nmolecules = # of non-zero IDs in molmap
// molmap[i] = index of molecule, skipping molecules not in group with -1
int
nmolecules
=
0
;
for
(
i
=
0
;
i
<
nlen
;
i
++
)
if
(
molmapall
[
i
])
molmap_H
[
i
]
=
nmolecules
++
;
else
molmap_H
[
i
]
=
-
1
;
memory
->
destroy
(
molmapall
);
// warn if any molecule has some atoms in group and some not in group
flag
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit
)
continue
;
if
(
molecule
[
i
]
<
idlo
||
molecule
[
i
]
>
idhi
)
continue
;
if
(
molmap_H
[
molecule
[
i
]
-
idlo
]
>=
0
)
flag
=
1
;
}
MPI_Allreduce
(
&
flag
,
&
flagall
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flagall
&&
comm
->
me
==
0
)
error
->
warning
(
FLERR
,
"One or more compute molecules has atoms not in group"
);
// if molmap simply stores 1 to Nmolecules, then free it
if
(
idlo
==
1
&&
idhi
==
nmolecules
&&
nlen
==
nmolecules
)
{
memory
->
destroy
(
molmap_H
);
molmap_H
=
NULL
;
}
return
nmolecules
;
}
void
FixLambdaHCalc
::
Print_Compensation_Density
(){
FILE
*
fp1
;
fp1
=
fopen
(
"Mean_Comp_Density.txt"
,
"w"
);
if
(
fp1
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open fix Mean_Comp_Density file %s"
,
"Mean_Comp_Density.txt"
);
error
->
one
(
FLERR
,
str
);
}
for
(
int
i
=
0
;
i
<
Density_Bin_Num
;
i
++
){
fprintf
(
fp1
,
"%d"
,
i
+
1
);
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
)
fprintf
(
fp1
,
"
\t
%.10f"
,
Mean_grad_Comp_Density_Conv_H
[
i
][
k
]);
fprintf
(
fp1
,
"
\n
"
);
}
fclose
(
fp1
);
/*
FILE *fp2;
char filename2[1000];
sprintf(filename2, "Mean_Comp_Density%4.1f_%d.txt", Comp_Density_Scaling_factor_H,Density_Update_Frequency);
fp2 = fopen(filename2,"a");
if (fp2 == NULL) {
char str[128];
sprintf(str,"Cannot open fix Mean_Comp_Density file %s","Mean_Comp_Density.txt");
error->one(FLERR,str);
}
printf("Density_Fluctuation = %f\n",Density_Fluctuation);
double x0;
if(Hybrid_Style==0)x0=x0lo;
else x0=0;
for(int i = 0;i < Density_Bin_Num; i++){
fprintf(fp2,"%d\t %.10f\t",i+1,i*Density_Bin_Size+x0);
for(int k = 0; k < atom->nmoltypesH;k++)
fprintf(fp2,"\t%.10f\t %.10f\t %.10f\t %.10f",Mean_Comp_Density_Conv_H[i][k],grad_Comp_Density_Conv_H[i][k],Mean_grad_Comp_Density_Conv_H[i][k]);
fprintf(fp2,"\n");
}
fclose(fp2);
FILE *fp3;
char filename3[1000];
sprintf(filename3, "Density_Fluctuation%4.1f_%d.txt", Comp_Density_Scaling_factor_H,Density_Update_Frequency);
fp3 = fopen(filename3,"a");
if (fp3 == NULL) {
char str[128];
sprintf(str,"Cannot open fix Density_Fluctuation file %s","Density_Fluctuation.txt");
error->one(FLERR,str);
}
int This_Step = update->ntimestep;
fprintf(fp3,"%d %.10f\n",This_Step,Density_Fluctuation);
fclose(fp3);
*/
}
void
FixLambdaHCalc
::
Clear_File_Compensation_Density
(){
FILE
*
fp1
;
fp1
=
fopen
(
"Mean_Comp_Density.txt"
,
"w"
);
if
(
fp1
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open fix Mean_Comp_Density file %s"
,
"Mean_Comp_Density.txt"
);
error
->
one
(
FLERR
,
str
);
}
fclose
(
fp1
);
/*
FILE *fp2;
char filename2[1000];
sprintf(filename2, "Mean_Comp_Density%4.1f_%d.txt", Comp_Density_Scaling_factor_H,Density_Update_Frequency);
fp2 = fopen(filename2,"w");
if (fp2 == NULL) {
char str[128];
sprintf(str,"Cannot open fix Mean_Comp_Density file %s","Mean_Comp_Density.txt");
error->one(FLERR,str);
}
fclose(fp2);
FILE *fp3;
char filename3[1000];
sprintf(filename3, "Density_Fluctuation%4.1f_%d.txt", Comp_Density_Scaling_factor_H,Density_Update_Frequency);
fp3 = fopen(filename3,"w");
if (fp3 == NULL) {
char str[128];
sprintf(str,"Cannot open fix Density_Fluctuation file %s","Density_Fluctuation.txt");
error->one(FLERR,str);
}
fclose(fp3);
*/
}
void
FixLambdaHCalc
::
Load_Compensation_Density
(){
if
(
me
==
0
){
FILE
*
fp1
;
char
str
[
128
];
fp1
=
fopen
(
"Mean_Comp_Density.txt"
,
"r"
);
if
(
fp1
==
NULL
)
{
sprintf
(
str
,
"Cannot open fix Mean_Comp_Density.txt file %s"
,
"Mean_Comp_Density.txt"
);
error
->
one
(
FLERR
,
str
);
}
int
i1
;
float
i2
;
while
(
!
feof
(
fp1
)){
fscanf
(
fp1
,
"%d"
,
&
i1
);
for
(
int
k
=
0
;
k
<
atom
->
nmoltypesH
;
k
++
){
fscanf
(
fp1
,
"
\t
%f"
,
&
i2
);
Mean_grad_Comp_Density_Conv_H
[
i1
-
1
][
k
]
=
(
double
)
i2
;
}
fscanf
(
fp1
,
"
\n
"
);
if
(
i1
>
Density_Bin_Num
){
sprintf
(
str
,
"Density bin number mismatches %d != %d"
,
Density_Bin_Num
,
i1
);
error
->
one
(
FLERR
,
str
);
}
}
fclose
(
fp1
);
if
(
me
==
0
){
if
(
screen
)
fprintf
(
screen
,
"
\n\n
Density componsation forces distributed successfully!
\n\n
"
);
if
(
logfile
)
fprintf
(
logfile
,
"
\n\n
Density componsation forces distributed successfully!
\n\n
"
);
}
}
MPI_Bcast
(
Mean_grad_Comp_Density_Conv_H
,
Density_Bin_Num
*
(
atom
->
nmoltypesH
+
1
),
MPI_DOUBLE
,
0
,
world
);
}
Event Timeline
Log In to Comment