本文整理汇总了C++中pr_rvecs函数的典型用法代码示例。如果您正苦于以下问题:C++ pr_rvecs函数的具体用法?C++ pr_rvecs怎么用?C++ pr_rvecs使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了pr_rvecs函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: list_xtc
void list_xtc(const char *fn)
{
t_fileio *xd;
int indent;
char buf[256];
rvec *x;
matrix box;
int nframe, natoms, step;
real prec, time;
gmx_bool bOK;
xd = open_xtc(fn, "r");
read_first_xtc(xd, &natoms, &step, &time, box, &x, &prec, &bOK);
nframe = 0;
do
{
sprintf(buf, "%s frame %d", fn, nframe);
indent = 0;
indent = pr_title(stdout, indent, buf);
pr_indent(stdout, indent);
fprintf(stdout, "natoms=%10d step=%10d time=%12.7e prec=%10g\n",
natoms, step, time, prec);
pr_rvecs(stdout, indent, "box", box, DIM);
pr_rvecs(stdout, indent, "x", x, natoms);
nframe++;
}
while (read_next_xtc(xd, natoms, &step, &time, box, x, &prec, &bOK));
if (!bOK)
{
fprintf(stderr, "\nWARNING: Incomplete frame at time %g\n", time);
}
sfree(x);
close_xtc(xd);
}
开发者ID:JehandadKhan,项目名称:gromacs,代码行数:35,代码来源:dump.c
示例2: check_box_unchanged
/* Check whether the box is unchanged */
static void check_box_unchanged(matrix f_box, matrix box, char fn[], warninp_t wi)
{
int i, ii;
gmx_bool bSame = TRUE;
char warn_buf[STRLEN];
for (i = 0; i < DIM; i++)
{
for (ii = 0; ii < DIM; ii++)
{
if (f_box[i][ii] != box[i][ii])
{
bSame = FALSE;
}
}
}
if (!bSame)
{
sprintf(warn_buf, "%s Box size in reference file %s differs from actual box size!",
RotStr, fn);
warning(wi, warn_buf);
pr_rvecs(stderr, 0, "Your box is:", box, 3);
pr_rvecs(stderr, 0, "Box in file:", f_box, 3);
}
}
开发者ID:alwanderer,项目名称:gromacs,代码行数:27,代码来源:readrot.c
示例3: sum_forces
static void sum_forces(int start,int end,rvec f[],rvec flr[])
{
int i;
if (gmx_debug_at) {
pr_rvecs(debug,0,"fsr",f+start,end-start);
pr_rvecs(debug,0,"flr",flr+start,end-start);
}
for(i=start; (i<end); i++)
rvec_inc(f[i],flr[i]);
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:11,代码来源:sim_util.c
示例4: list_trn
static void list_trn(char *fn)
{
int fpread,fpwrite,nframe,indent;
char buf[256];
rvec *x,*v,*f;
matrix box;
t_trnheader trn;
bool bOK;
fpread = open_trn(fn,"r");
fpwrite = open_tpx(NULL,"w");
gmx_fio_setdebug(fpwrite,TRUE);
nframe = 0;
while (fread_trnheader(fpread,&trn,&bOK)) {
snew(x,trn.natoms);
snew(v,trn.natoms);
snew(f,trn.natoms);
if (fread_htrn(fpread,&trn,
trn.box_size ? box : NULL,
trn.x_size ? x : NULL,
trn.v_size ? v : NULL,
trn.f_size ? f : NULL)) {
sprintf(buf,"%s frame %d",fn,nframe);
indent=0;
indent=pr_title(stdout,indent,buf);
pr_indent(stdout,indent);
fprintf(stdout,"natoms=%10d step=%10d time=%12.7e lambda=%10g\n",
trn.natoms,trn.step,trn.t,trn.lambda);
if (trn.box_size)
pr_rvecs(stdout,indent,"box",box,DIM);
if (trn.x_size)
pr_rvecs(stdout,indent,"x",x,trn.natoms);
if (trn.v_size)
pr_rvecs(stdout,indent,"v",v,trn.natoms);
if (trn.f_size)
pr_rvecs(stdout,indent,"f",f,trn.natoms);
}
else
fprintf(stderr,"\nWARNING: Incomplete frame: nr %d, t=%g\n",
nframe,trn.t);
sfree(x);
sfree(v);
sfree(f);
nframe++;
}
if (!bOK)
fprintf(stderr,"\nWARNING: Incomplete frame header: nr %d, t=%g\n",
nframe,trn.t);
close_tpx(fpwrite);
close_trn(fpread);
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:53,代码来源:gmxdump.c
示例5: dump_it_all
static void dump_it_all(FILE *fp,char *title,
int natoms,rvec x[],rvec xp[],rvec v[],rvec f[])
{
#ifdef DEBUG
if (fp) {
fprintf(fp,"%s\n",title);
pr_rvecs(fp,0,"x",x,natoms);
pr_rvecs(fp,0,"xp",xp,natoms);
pr_rvecs(fp,0,"v",v,natoms);
pr_rvecs(fp,0,"f",f,natoms);
}
#endif
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:13,代码来源:update.c
示例6: correct_box_elem
static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
{
int shift, maxshift = 10;
shift = 0;
/* correct elem d of vector v with vector d */
while (box[v][d] > BOX_MARGIN_CORRECT*0.5*box[d][d])
{
if (fplog)
{
fprintf(fplog, "Step %d: correcting invalid box:\n", step);
pr_rvecs(fplog, 0, "old box", box, DIM);
}
rvec_dec(box[v], box[d]);
shift--;
if (fplog)
{
pr_rvecs(fplog, 0, "new box", box, DIM);
}
if (shift <= -maxshift)
{
gmx_fatal(FARGS,
"Box was shifted at least %d times. Please see log-file.",
maxshift);
}
}
while (box[v][d] < -BOX_MARGIN_CORRECT*0.5*box[d][d])
{
if (fplog)
{
fprintf(fplog, "Step %d: correcting invalid box:\n", step);
pr_rvecs(fplog, 0, "old box", box, DIM);
}
rvec_inc(box[v], box[d]);
shift++;
if (fplog)
{
pr_rvecs(fplog, 0, "new box", box, DIM);
}
if (shift >= maxshift)
{
gmx_fatal(FARGS,
"Box was shifted at least %d times. Please see log-file.",
maxshift);
}
}
return shift;
}
开发者ID:yhalcyon,项目名称:Gromacs,代码行数:50,代码来源:pbc.c
示例7: calc_force
void calc_force(int natom,rvec f[],rvec fff[])
{
int i,j,m;
int jindex[] = { 0, 5, 10};
rvec dx,df;
real msf1,msf2;
for(j=0; (j<2); j++) {
clear_rvec(fff[j]);
for(i=jindex[j]; (i<jindex[j+1]); i++) {
for(m=0; (m<DIM); m++) {
fff[j][m] += f[i][m];
}
}
}
msf1 = iprod(fff[0],fff[0]);
msf2 = iprod(fff[1],fff[1]);
if (debug) {
pr_rvecs(debug,0,"force",f,natom);
fprintf(debug,"FMOL: %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n",
fff[0][XX],fff[0][YY],fff[0][ZZ],fff[1][XX],fff[1][YY],fff[1][ZZ]);
fprintf(debug,"RMSF: %10.3e %10.3e\n",msf1,msf2);
}
}
开发者ID:ElsevierSoftwareX,项目名称:SOFTX-D-15-00003,代码行数:26,代码来源:calcfdev.c
示例8: calc_virial
static void calc_virial(FILE *fplog,int start,int homenr,rvec x[],rvec f[],
tensor vir_part,t_graph *graph,matrix box,
t_nrnb *nrnb,const t_forcerec *fr,int ePBC)
{
int i,j;
tensor virtest;
/* The short-range virial from surrounding boxes */
clear_mat(vir_part);
calc_vir(fplog,SHIFTS,fr->shift_vec,fr->fshift,vir_part,ePBC==epbcSCREW,box);
inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
/* Calculate partial virial, for local atoms only, based on short range.
* Total virial is computed in global_stat, called from do_md
*/
f_calc_vir(fplog,start,start+homenr,x,f,vir_part,graph,box);
inc_nrnb(nrnb,eNR_VIRIAL,homenr);
/* Add position restraint contribution */
for(i=0; i<DIM; i++) {
vir_part[i][i] += fr->vir_diag_posres[i];
}
/* Add wall contribution */
for(i=0; i<DIM; i++) {
vir_part[i][ZZ] += fr->vir_wall_z[i];
}
if (debug)
pr_rvecs(debug,0,"vir_part",vir_part,DIM);
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:31,代码来源:sim_util.c
示例9: pr_matrix
static void pr_matrix(FILE *fp,int indent,const char *title,rvec *m,
gmx_bool bMDPformat)
{
if (bMDPformat)
fprintf(fp,"%-10s = %g %g %g %g %g %g\n",title,
m[XX][XX],m[YY][YY],m[ZZ][ZZ],m[XX][YY],m[XX][ZZ],m[YY][ZZ]);
else
pr_rvecs(fp,indent,title,m,DIM);
}
开发者ID:alexholehouse,项目名称:gromacs,代码行数:9,代码来源:txtdump.c
示例10: init_parallel
void init_parallel(FILE *log,char *tpxfile,t_commrec *cr,
t_inputrec *inputrec,gmx_mtop_t *mtop,
t_state *state,
int list)
{
int step;
real t;
char buf[256];
if (MASTER(cr)) {
init_inputrec(inputrec);
read_tpx_state(tpxfile,&step,&t,inputrec,state,NULL,mtop);
/* When we will be doing domain decomposition with separate PME nodes
* the rng entries will be too large, we correct for this later.
*/
set_state_entries(state,inputrec,cr->nnodes);
}
bcast_ir_mtop(cr,inputrec,mtop);
if (inputrec->eI == eiBD || EI_SD(inputrec->eI)) {
/* Make sure the random seeds are different on each node */
inputrec->ld_seed += cr->nodeid;
}
/* Printing */
if (list!=0 && log!=NULL)
{
if (list&LIST_INPUTREC)
pr_inputrec(log,0,"parameters of the run",inputrec,FALSE);
if (list&LIST_X)
pr_rvecs(log,0,"box",state->box,DIM);
if (list&LIST_X)
pr_rvecs(log,0,"box_rel",state->box_rel,DIM);
if (list&LIST_V)
pr_rvecs(log,0,"boxv",state->boxv,DIM);
if (list&LIST_X)
pr_rvecs(log,0,int_title("x",0,buf,255),state->x,state->natoms);
if (list&LIST_V)
pr_rvecs(log,0,int_title("v",0,buf,255),state->v,state->natoms);
if (list&LIST_TOP)
pr_mtop(log,0,int_title("topology",cr->nodeid,buf,255),mtop,TRUE);
fflush(log);
}
}
开发者ID:BioinformaticsArchive,项目名称:GromPy,代码行数:44,代码来源:init.c
示例11: pr_molblock
static void pr_molblock(FILE *fp,int indent,const char *title,
gmx_molblock_t *molb,int n,
gmx_moltype_t *molt,
gmx_bool bShowNumbers)
{
indent = pr_title_n(fp,indent,title,n);
(void) pr_indent(fp,indent);
(void) fprintf(fp,"%-20s = %d \"%s\"\n",
"moltype",molb->type,*(molt[molb->type].name));
pr_int(fp,indent,"#molecules",molb->nmol);
pr_int(fp,indent,"#atoms_mol",molb->natoms_mol);
pr_int(fp,indent,"#posres_xA",molb->nposres_xA);
if (molb->nposres_xA > 0) {
pr_rvecs(fp,indent,"posres_xA",molb->posres_xA,molb->nposres_xA);
}
pr_int(fp,indent,"#posres_xB",molb->nposres_xB);
if (molb->nposres_xB > 0) {
pr_rvecs(fp,indent,"posres_xB",molb->posres_xB,molb->nposres_xB);
}
}
开发者ID:alexholehouse,项目名称:gromacs,代码行数:20,代码来源:txtdump.c
示例12: check_pbc
static void check_pbc(FILE *fp,rvec x[],int shell)
{
int m,now;
now = shell-4;
for(m=0; (m<DIM); m++)
if (fabs(x[shell][m]-x[now][m]) > 0.3) {
pr_rvecs(fp,0,"SHELL-X",x+now,5);
break;
}
}
开发者ID:martinhoefling,项目名称:gromacs,代码行数:11,代码来源:shellfc.c
示例13: list_xtc
void list_xtc(char *fn, bool bXVG)
{
int xd,indent;
char buf[256];
rvec *x;
matrix box;
int nframe,natoms,step;
real prec,time;
bool bOK;
xd = open_xtc(fn,"r");
read_first_xtc(xd,&natoms,&step,&time,box,&x,&prec,&bOK);
nframe=0;
do {
if (bXVG) {
int i,d;
fprintf(stdout,"%g",time);
for(i=0; i<natoms; i++)
for(d=0; d<DIM; d++)
fprintf(stdout," %g",x[i][d]);
fprintf(stdout,"\n");
} else {
sprintf(buf,"%s frame %d",fn,nframe);
indent=0;
indent=pr_title(stdout,indent,buf);
pr_indent(stdout,indent);
fprintf(stdout,"natoms=%10d step=%10d time=%12.7e prec=%10g\n",
natoms,step,time,prec);
pr_rvecs(stdout,indent,"box",box,DIM);
pr_rvecs(stdout,indent,"x",x,natoms);
}
nframe++;
} while (read_next_xtc(xd,natoms,&step,&time,box,x,&prec,&bOK));
if (!bOK)
fprintf(stderr,"\nWARNING: Incomplete frame at time %g\n",time);
close_xtc(xd);
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:39,代码来源:gmxdump.c
示例14: calc_gyro
real calc_gyro(rvec x[], int gnx, atom_id index[], t_atom atom[], real tm,
rvec gvec, rvec d, gmx_bool bQ, gmx_bool bRot, gmx_bool bMOI, matrix trans)
{
int i, ii, m;
real gyro, dx2, m0, Itot;
rvec comp;
if (bRot)
{
principal_comp(gnx, index, atom, x, trans, d);
Itot = norm(d);
if (bMOI)
{
return Itot;
}
for (m = 0; (m < DIM); m++)
{
d[m] = std::sqrt(d[m]/tm);
}
#ifdef DEBUG
pr_rvecs(stderr, 0, "trans", trans, DIM);
#endif
/* rotate_atoms(gnx,index,x,trans); */
}
clear_rvec(comp);
for (i = 0; (i < gnx); i++)
{
ii = index[i];
if (bQ)
{
m0 = std::abs(atom[ii].q);
}
else
{
m0 = atom[ii].m;
}
for (m = 0; (m < DIM); m++)
{
dx2 = x[ii][m]*x[ii][m];
comp[m] += dx2*m0;
}
}
gyro = comp[XX]+comp[YY]+comp[ZZ];
for (m = 0; (m < DIM); m++)
{
gvec[m] = std::sqrt((gyro-comp[m])/tm);
}
return std::sqrt(gyro/tm);
}
开发者ID:tanigawa,项目名称:gromacs,代码行数:51,代码来源:gmx_gyrate.cpp
示例15: calc_virial
void calc_virial(FILE *log,int start,int homenr,rvec x[],rvec f[],
tensor vir_part,tensor pme_vir,
t_graph *graph,matrix box,
t_nrnb *nrnb,t_forcerec *fr,bool bTweak)
{
int i,j;
tensor virtest;
/* Now it is time for the short range virial. At this timepoint vir_part
* already contains the virial from surrounding boxes.
* Calculate partial virial, for local atoms only, based on short range.
* Total virial is computed in global_stat, called from do_md
*/
f_calc_vir(log,start,start+homenr,x,f,vir_part,graph,box);
inc_nrnb(nrnb,eNR_VIRIAL,homenr);
/* Add up the long range forces if necessary */
/* if (!bTweak) {
sum_lrforces(f,fr,start,homenr);
}*/
/* Add up virial if necessary */
if (EEL_LR(fr->eeltype) && (fr->eeltype != eelPPPM)) {
if (debug && bTweak) {
clear_mat(virtest);
f_calc_vir(log,start,start+homenr,x,fr->f_pme,virtest,graph,box);
pr_rvecs(debug,0,"virtest",virtest,DIM);
pr_rvecs(debug,0,"pme_vir",pme_vir,DIM);
}
/* PPPM virial sucks */
if (!bTweak)
for(i=0; (i<DIM); i++)
for(j=0; (j<DIM); j++)
vir_part[i][j]+=pme_vir[i][j];
}
if (debug)
pr_rvecs(debug,0,"vir_part",vir_part,DIM);
}
开发者ID:Chadi-akel,项目名称:cere,代码行数:38,代码来源:sim_util.c
示例16: calc_angles_dihs
void calc_angles_dihs(t_params *ang, t_params *dih, rvec x[], gmx_bool bPBC,
matrix box)
{
int i, ai, aj, ak, al, t1, t2, t3;
rvec r_ij, r_kj, r_kl, m, n;
real sign, th, costh, ph;
t_pbc pbc;
if (bPBC)
{
set_pbc(&pbc, epbcXYZ, box);
}
if (debug)
{
pr_rvecs(debug, 0, "X2TOP", box, DIM);
}
for (i = 0; (i < ang->nr); i++)
{
ai = ang->param[i].AI;
aj = ang->param[i].AJ;
ak = ang->param[i].AK;
th = RAD2DEG*bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : NULL,
r_ij, r_kj, &costh, &t1, &t2);
if (debug)
{
fprintf(debug, "X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
ai, aj, ak, norm(r_ij), norm(r_kj), th);
}
ang->param[i].C0 = th;
}
for (i = 0; (i < dih->nr); i++)
{
ai = dih->param[i].AI;
aj = dih->param[i].AJ;
ak = dih->param[i].AK;
al = dih->param[i].AL;
ph = RAD2DEG*dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : NULL,
r_ij, r_kj, r_kl, m, n, &sign, &t1, &t2, &t3);
if (debug)
{
fprintf(debug, "X2TOP: ai=%3d aj=%3d ak=%3d al=%3d r_ij=%8.3f r_kj=%8.3f r_kl=%8.3f ph=%8.3f\n",
ai, aj, ak, al, norm(r_ij), norm(r_kj), norm(r_kl), ph);
}
dih->param[i].C0 = ph;
}
}
开发者ID:exianshine,项目名称:gromacs,代码行数:46,代码来源:g_x2top.c
示例17: pr_rotgrp
static void pr_rotgrp(FILE *fp,int indent,int g,t_rotgrp *rotg)
{
pr_indent(fp,indent);
fprintf(fp,"rotation_group %d:\n",g);
indent += 2;
PS("type",EROTGEOM(rotg->eType));
PS("massw",BOOL(rotg->bMassW));
pr_ivec_block(fp,indent,"atom",rotg->ind,rotg->nat,TRUE);
pr_rvecs(fp,indent,"x_ref",rotg->x_ref,rotg->nat);
pr_rvec(fp,indent,"vec",rotg->vec,DIM,TRUE);
pr_rvec(fp,indent,"pivot",rotg->pivot,DIM,TRUE);
PR("rate",rotg->rate);
PR("k",rotg->k);
PR("slab_dist",rotg->slab_dist);
PR("min_gaussian",rotg->min_gaussian);
PR("epsilon",rotg->eps);
PS("fit_method",EROTFIT(rotg->eFittype));
PI("potfitangle_nstep",rotg->PotAngle_nstep);
PR("potfitangle_step",rotg->PotAngle_step);
}
开发者ID:alexholehouse,项目名称:gromacs,代码行数:20,代码来源:txtdump.c
示例18: dump_pbc
void dump_pbc(FILE *fp, t_pbc *pbc)
{
rvec sum_box;
fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
pr_rvecs(fp, 0, "box", pbc->box, DIM);
pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
if (pbc->ntric_vec > 0)
{
pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
}
}
开发者ID:wangxubo0201,项目名称:gromacs,代码行数:19,代码来源:pbc.cpp
示例19: list_trr
static void list_trr(const char *fn)
{
t_fileio *fpread;
int nframe, indent;
char buf[256];
rvec *x, *v, *f;
matrix box;
gmx_trr_header_t trrheader;
gmx_bool bOK;
fpread = gmx_trr_open(fn, "r");
nframe = 0;
while (gmx_trr_read_frame_header(fpread, &trrheader, &bOK))
{
snew(x, trrheader.natoms);
snew(v, trrheader.natoms);
snew(f, trrheader.natoms);
if (gmx_trr_read_frame_data(fpread, &trrheader,
trrheader.box_size ? box : NULL,
trrheader.x_size ? x : NULL,
trrheader.v_size ? v : NULL,
trrheader.f_size ? f : NULL))
{
sprintf(buf, "%s frame %d", fn, nframe);
indent = 0;
indent = pr_title(stdout, indent, buf);
pr_indent(stdout, indent);
fprintf(stdout, "natoms=%10d step=%10d time=%12.7e lambda=%10g\n",
trrheader.natoms, trrheader.step, trrheader.t, trrheader.lambda);
if (trrheader.box_size)
{
pr_rvecs(stdout, indent, "box", box, DIM);
}
if (trrheader.x_size)
{
pr_rvecs(stdout, indent, "x", x, trrheader.natoms);
}
if (trrheader.v_size)
{
pr_rvecs(stdout, indent, "v", v, trrheader.natoms);
}
if (trrheader.f_size)
{
pr_rvecs(stdout, indent, "f", f, trrheader.natoms);
}
}
else
{
fprintf(stderr, "\nWARNING: Incomplete frame: nr %d, t=%g\n",
nframe, trrheader.t);
}
sfree(x);
sfree(v);
sfree(f);
nframe++;
}
if (!bOK)
{
fprintf(stderr, "\nWARNING: Incomplete frame header: nr %d, t=%g\n",
nframe, trrheader.t);
}
gmx_trr_close(fpread);
}
开发者ID:carryer123,项目名称:gromacs,代码行数:65,代码来源:dump.cpp
示例20: check_cm_grp
//.........这里部分代码省略.........
if (vcm->mode != ecmNO)
{
for (g = 0; (g < vcm->nr); g++)
{
tm = vcm->group_mass[g];
if (tm != 0)
{
tm_1 = 1.0/tm;
svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
}
/* Else it's zero anyway! */
}
if (vcm->mode == ecmANGULAR)
{
for (g = 0; (g < vcm->nr); g++)
{
tm = vcm->group_mass[g];
if (tm != 0)
{
tm_1 = 1.0/tm;
/* Compute center of mass for this group */
for (m = 0; (m < DIM); m++)
{
vcm->group_x[g][m] *= tm_1;
}
/* Subtract the center of mass contribution to the
* angular momentum
*/
cprod(vcm->group_x[g], vcm->group_v[g], jcm);
for (m = 0; (m < DIM); m++)
{
vcm->group_j[g][m] -= tm*jcm[m];
}
/* Subtract the center of mass contribution from the inertia
* tensor (this is not as trivial as it seems, but due to
* some cancellation we can still do it, even in parallel).
*/
clear_mat(Icm);
update_tensor(vcm->group_x[g], tm, Icm);
m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
/* Compute angular velocity, using matrix operation
* Since J = I w
* we have
* w = I^-1 J
*/
get_minv(vcm->group_i[g], Icm);
mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
}
/* Else it's zero anyway! */
}
}
}
for (g = 0; (g < vcm->nr); g++)
{
ekcm = 0;
if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
{
for (m = 0; m < vcm->ndim; m++)
{
ekcm += sqr(vcm->group_v[g][m]);
}
ekcm *= 0.5*vcm->group_mass[g];
Temp_cm = 2*ekcm/vcm->group_ndf[g];
if ((Temp_cm > Temp_Max) && fp)
{
fprintf(fp, "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
vcm->group_name[g], vcm->group_v[g][XX],
vcm->group_v[g][YY], vcm->group_v[g][ZZ], Temp_cm);
}
if (vcm->mode == ecmANGULAR)
{
ekrot = 0.5*iprod(vcm->group_j[g], vcm->group_w[g]);
if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI))
{
/* if we have an integrator that may not conserve momenta, skip */
tm = vcm->group_mass[g];
fprintf(fp, "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
vcm->group_name[g], tm, ekrot, det(vcm->group_i[g]));
fprintf(fp, " COM: %12.5f %12.5f %12.5f\n",
vcm->group_x[g][XX], vcm->group_x[g][YY], vcm->group_x[g][ZZ]);
fprintf(fp, " P: %12.5f %12.5f %12.5f\n",
vcm->group_p[g][XX], vcm->group_p[g][YY], vcm->group_p[g][ZZ]);
fprintf(fp, " V: %12.5f %12.5f %12.5f\n",
vcm->group_v[g][XX], vcm->group_v[g][YY], vcm->group_v[g][ZZ]);
fprintf(fp, " J: %12.5f %12.5f %12.5f\n",
vcm->group_j[g][XX], vcm->group_j[g][YY], vcm->group_j[g][ZZ]);
fprintf(fp, " w: %12.5f %12.5f %12.5f\n",
vcm->group_w[g][XX], vcm->group_w[g][YY], vcm->group_w[g][ZZ]);
pr_rvecs(fp, 0, "Inertia tensor", vcm->group_i[g], DIM);
}
}
}
}
}
开发者ID:rmcgibbo,项目名称:gromacs,代码行数:101,代码来源:vcm.cpp
注:本文中的pr_rvecs函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论