Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71663688
fix_spring_chunk.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
Fri, Jul 12, 12:24
Size
6 KB
Mime Type
text/x-c
Expires
Sun, Jul 14, 12:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18967572
Attached To
rLAMMPS lammps
fix_spring_chunk.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.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "fix_spring_chunk.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "respa.h"
#include "domain.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "compute_com_chunk.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
#define SMALL 1.0e-10
/* ---------------------------------------------------------------------- */
FixSpringChunk
::
FixSpringChunk
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
)
{
if
(
narg
!=
6
)
error
->
all
(
FLERR
,
"Illegal fix spring/chunk command"
);
scalar_flag
=
1
;
global_freq
=
1
;
extscalar
=
1
;
k_spring
=
force
->
numeric
(
FLERR
,
arg
[
3
]);
int
n
=
strlen
(
arg
[
4
])
+
1
;
idchunk
=
new
char
[
n
];
strcpy
(
idchunk
,
arg
[
4
]);
n
=
strlen
(
arg
[
5
])
+
1
;
idcom
=
new
char
[
n
];
strcpy
(
idcom
,
arg
[
5
]);
esprings
=
0.0
;
com0
=
fcom
=
NULL
;
nchunk
=
0
;
}
/* ---------------------------------------------------------------------- */
FixSpringChunk
::~
FixSpringChunk
()
{
memory
->
destroy
(
com0
);
memory
->
destroy
(
fcom
);
// decrement lock counter in compute chunk/atom, it if still exists
int
icompute
=
modify
->
find_compute
(
idchunk
);
if
(
icompute
>=
0
)
{
cchunk
=
(
ComputeChunkAtom
*
)
modify
->
compute
[
icompute
];
cchunk
->
unlock
(
this
);
cchunk
->
lockcount
--
;
}
delete
[]
idchunk
;
delete
[]
idcom
;
}
/* ---------------------------------------------------------------------- */
int
FixSpringChunk
::
setmask
()
{
int
mask
=
0
;
mask
|=
POST_FORCE
;
mask
|=
THERMO_ENERGY
;
mask
|=
POST_FORCE_RESPA
;
mask
|=
MIN_POST_FORCE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixSpringChunk
::
init
()
{
// current indices for idchunk and idcom
int
icompute
=
modify
->
find_compute
(
idchunk
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Chunk/atom compute does not exist for fix spring/chunk"
);
cchunk
=
(
ComputeChunkAtom
*
)
modify
->
compute
[
icompute
];
if
(
strcmp
(
cchunk
->
style
,
"chunk/atom"
)
!=
0
)
error
->
all
(
FLERR
,
"Fix spring/chunk does not use chunk/atom compute"
);
icompute
=
modify
->
find_compute
(
idcom
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Com/chunk compute does not exist for fix spring/chunk"
);
ccom
=
(
ComputeCOMChunk
*
)
modify
->
compute
[
icompute
];
if
(
strcmp
(
ccom
->
style
,
"com/chunk"
)
!=
0
)
error
->
all
(
FLERR
,
"Fix spring/chunk does not use com/chunk compute"
);
// check that idchunk is consistent with ccom->idchunk
if
(
strcmp
(
idchunk
,
ccom
->
idchunk
)
!=
0
)
error
->
all
(
FLERR
,
"Fix spring chunk chunkID not same as comID chunkID"
);
if
(
strstr
(
update
->
integrate_style
,
"respa"
))
nlevels_respa
=
((
Respa
*
)
update
->
integrate
)
->
nlevels
;
}
/* ---------------------------------------------------------------------- */
void
FixSpringChunk
::
setup
(
int
vflag
)
{
if
(
strstr
(
update
->
integrate_style
,
"verlet"
))
post_force
(
vflag
);
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
FixSpringChunk
::
min_setup
(
int
vflag
)
{
post_force
(
vflag
);
}
/* ---------------------------------------------------------------------- */
void
FixSpringChunk
::
post_force
(
int
vflag
)
{
int
i
,
m
;
double
dx
,
dy
,
dz
,
r
;
// check if first time cchunk will be queried via ccom
// if so, lock idchunk for as long as this fix is in place
// will be unlocked in destructor
// necessary b/c this fix stores original COM
if
(
com0
==
NULL
)
cchunk
->
lock
(
this
,
update
->
ntimestep
,
-
1
);
// calculate current centers of mass for each chunk
// extract pointers from idchunk and idcom
ccom
->
compute_array
();
nchunk
=
cchunk
->
nchunk
;
int
*
ichunk
=
cchunk
->
ichunk
;
double
*
masstotal
=
ccom
->
masstotal
;
double
**
com
=
ccom
->
array
;
// check if first time cchunk was queried via ccom
// if so, allocate com0,fcom and store initial COM
if
(
com0
==
NULL
)
{
memory
->
create
(
com0
,
nchunk
,
3
,
"spring/chunk:com0"
);
memory
->
create
(
fcom
,
nchunk
,
3
,
"spring/chunk:fcom"
);
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
com0
[
m
][
0
]
=
com
[
m
][
0
];
com0
[
m
][
1
]
=
com
[
m
][
1
];
com0
[
m
][
2
]
=
com
[
m
][
2
];
}
}
// calculate fcom = force on each COM, divided by masstotal
esprings
=
0.0
;
for
(
m
=
0
;
m
<
nchunk
;
m
++
)
{
dx
=
com
[
m
][
0
]
-
com0
[
m
][
0
];
dy
=
com
[
m
][
1
]
-
com0
[
m
][
1
];
dz
=
com
[
m
][
2
]
-
com0
[
m
][
2
];
r
=
sqrt
(
dx
*
dx
+
dy
*
dy
+
dz
*
dz
);
r
=
MAX
(
r
,
SMALL
);
if
(
masstotal
[
m
])
{
fcom
[
m
][
0
]
=
k_spring
*
dx
/
r
/
masstotal
[
m
];
fcom
[
m
][
1
]
=
k_spring
*
dy
/
r
/
masstotal
[
m
];
fcom
[
m
][
2
]
=
k_spring
*
dz
/
r
/
masstotal
[
m
];
esprings
+=
0.5
*
k_spring
*
r
*
r
;
}
}
// apply restoring force to atoms in each chunk
double
**
f
=
atom
->
f
;
int
*
type
=
atom
->
type
;
double
*
mass
=
atom
->
mass
;
double
*
rmass
=
atom
->
rmass
;
int
nlocal
=
atom
->
nlocal
;
double
massone
;
if
(
rmass
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
m
=
ichunk
[
i
]
-
1
;
if
(
m
<
0
)
continue
;
massone
=
rmass
[
i
];
f
[
i
][
0
]
-=
fcom
[
m
][
0
]
*
massone
;
f
[
i
][
1
]
-=
fcom
[
m
][
1
]
*
massone
;
f
[
i
][
2
]
-=
fcom
[
m
][
2
]
*
massone
;
}
}
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
m
=
ichunk
[
i
]
-
1
;
if
(
m
<
0
)
continue
;
massone
=
mass
[
type
[
i
]];
f
[
i
][
0
]
-=
fcom
[
m
][
0
]
*
massone
;
f
[
i
][
1
]
-=
fcom
[
m
][
1
]
*
massone
;
f
[
i
][
2
]
-=
fcom
[
m
][
2
]
*
massone
;
}
}
}
/* ---------------------------------------------------------------------- */
void
FixSpringChunk
::
post_force_respa
(
int
vflag
,
int
ilevel
,
int
iloop
)
{
if
(
ilevel
==
nlevels_respa
-
1
)
post_force
(
vflag
);
}
/* ---------------------------------------------------------------------- */
void
FixSpringChunk
::
min_post_force
(
int
vflag
)
{
post_force
(
vflag
);
}
/* ----------------------------------------------------------------------
energy of stretched spring
------------------------------------------------------------------------- */
double
FixSpringChunk
::
compute_scalar
()
{
return
esprings
;
}
Event Timeline
Log In to Comment