本文整理汇总了C++中copy_rvec函数的典型用法代码示例。如果您正苦于以下问题:C++ copy_rvec函数的具体用法?C++ copy_rvec怎么用?C++ copy_rvec使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了copy_rvec函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: make_cyl_refgrps
static void make_cyl_refgrps(t_commrec *cr, struct pull_t *pull, t_mdatoms *md,
t_pbc *pbc, double t, rvec *x)
{
/* The size and stride per coord for the reduction buffer */
const int stride = 9;
int c, i, ii, m, start, end;
rvec g_x, dx, dir;
double inv_cyl_r2;
pull_comm_t *comm;
gmx_ga2la_t *ga2la = NULL;
comm = &pull->comm;
if (comm->dbuf_cyl == NULL)
{
snew(comm->dbuf_cyl, pull->ncoord*stride);
}
if (cr && DOMAINDECOMP(cr))
{
ga2la = cr->dd->ga2la;
}
start = 0;
end = md->homenr;
inv_cyl_r2 = 1.0/gmx::square(pull->params.cylinder_r);
/* loop over all groups to make a reference group for each*/
for (c = 0; c < pull->ncoord; c++)
{
pull_coord_work_t *pcrd;
double sum_a, wmass, wwmass;
dvec radf_fac0, radf_fac1;
pcrd = &pull->coord[c];
sum_a = 0;
wmass = 0;
wwmass = 0;
clear_dvec(radf_fac0);
clear_dvec(radf_fac1);
if (pcrd->params.eGeom == epullgCYL)
{
pull_group_work_t *pref, *pgrp, *pdyna;
/* pref will be the same group for all pull coordinates */
pref = &pull->group[pcrd->params.group[0]];
pgrp = &pull->group[pcrd->params.group[1]];
pdyna = &pull->dyna[c];
copy_rvec(pcrd->vec, dir);
pdyna->nat_loc = 0;
/* We calculate distances with respect to the reference location
* of this cylinder group (g_x), which we already have now since
* we reduced the other group COM over the ranks. This resolves
* any PBC issues and we don't need to use a PBC-atom here.
*/
if (pcrd->params.rate != 0)
{
/* With rate=0, value_ref is set initially */
pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
}
for (m = 0; m < DIM; m++)
{
g_x[m] = pgrp->x[m] - pcrd->vec[m]*pcrd->value_ref;
}
/* loop over all atoms in the main ref group */
for (i = 0; i < pref->params.nat; i++)
{
ii = pref->params.ind[i];
if (ga2la)
{
if (!ga2la_get_home(ga2la, pref->params.ind[i], &ii))
{
ii = -1;
}
}
if (ii >= start && ii < end)
{
double dr2, dr2_rel, inp;
dvec dr;
pbc_dx_aiuc(pbc, x[ii], g_x, dx);
inp = iprod(dir, dx);
dr2 = 0;
for (m = 0; m < DIM; m++)
{
/* Determine the radial components */
dr[m] = dx[m] - inp*dir[m];
dr2 += dr[m]*dr[m];
}
dr2_rel = dr2*inv_cyl_r2;
if (dr2_rel < 1)
{
double mass, weight, dweight_r;
dvec mdw;
//.........这里部分代码省略.........
开发者ID:MelroLeandro,项目名称:gromacs,代码行数:101,代码来源:pullutil.cpp
示例2: do_tpi
//.........这里部分代码省略.........
atoms2md(top_global, inputrec, 0, NULL, top_global->natoms, mdatoms);
update_mdatoms(mdatoms, inputrec->fepvals->init_lambda);
snew(enerd, 1);
init_enerdata(groups->grps[egcENER].nr, inputrec->fepvals->n_lambda, enerd);
snew(f, top_global->natoms);
/* Print to log file */
walltime_accounting_start(walltime_accounting);
wallcycle_start(wcycle, ewcRUN);
print_start(fplog, cr, walltime_accounting, "Test Particle Insertion");
/* The last charge group is the group to be inserted */
cg_tp = top->cgs.nr - 1;
a_tp0 = top->cgs.index[cg_tp];
a_tp1 = top->cgs.index[cg_tp+1];
if (debug)
{
fprintf(debug, "TPI cg %d, atoms %d-%d\n", cg_tp, a_tp0, a_tp1);
}
if (a_tp1 - a_tp0 > 1 &&
(inputrec->rlist < inputrec->rcoulomb ||
inputrec->rlist < inputrec->rvdw))
{
gmx_fatal(FARGS, "Can not do TPI for multi-atom molecule with a twin-range cut-off");
}
snew(x_mol, a_tp1-a_tp0);
bDispCorr = (inputrec->eDispCorr != edispcNO);
bCharge = FALSE;
for (i = a_tp0; i < a_tp1; i++)
{
/* Copy the coordinates of the molecule to be insterted */
copy_rvec(state->x[i], x_mol[i-a_tp0]);
/* Check if we need to print electrostatic energies */
bCharge |= (mdatoms->chargeA[i] != 0 ||
(mdatoms->chargeB && mdatoms->chargeB[i] != 0));
}
bRFExcl = (bCharge && EEL_RF(fr->eeltype) && fr->eeltype != eelRF_NEC);
calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), state->x, fr->cg_cm);
if (bCavity)
{
if (norm(fr->cg_cm[cg_tp]) > 0.5*inputrec->rlist && fplog)
{
fprintf(fplog, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
fprintf(stderr, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
}
}
else
{
/* Center the molecule to be inserted at zero */
for (i = 0; i < a_tp1-a_tp0; i++)
{
rvec_dec(x_mol[i], fr->cg_cm[cg_tp]);
}
}
if (fplog)
{
fprintf(fplog, "\nWill insert %d atoms %s partial charges\n",
a_tp1-a_tp0, bCharge ? "with" : "without");
fprintf(fplog, "\nWill insert %d times in each frame of %s\n",
(int)nsteps, opt2fn("-rerun", nfile, fnm));
}
开发者ID:daniellandau,项目名称:gromacs,代码行数:67,代码来源:tpi.c
示例3: write_trxframe_indexed
int write_trxframe_indexed(t_trxstatus *status, const t_trxframe *fr, int nind,
const int *ind, gmx_conect gc)
{
char title[STRLEN];
rvec *xout = NULL, *vout = NULL, *fout = NULL;
int i, ftp = -1;
real prec;
if (fr->bPrec)
{
prec = fr->prec;
}
else
{
prec = 1000.0;
}
if (status->tng)
{
ftp = efTNG;
}
else if (status->fio)
{
ftp = gmx_fio_getftp(status->fio);
}
else
{
gmx_incons("No input file available");
}
switch (ftp)
{
case efTRR:
case efTNG:
break;
default:
if (!fr->bX)
{
gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
ftp2ext(ftp));
}
break;
}
switch (ftp)
{
case efTRR:
case efTNG:
if (fr->bV)
{
snew(vout, nind);
for (i = 0; i < nind; i++)
{
copy_rvec(fr->v[ind[i]], vout[i]);
}
}
if (fr->bF)
{
snew(fout, nind);
for (i = 0; i < nind; i++)
{
copy_rvec(fr->f[ind[i]], fout[i]);
}
}
/* no break */
case efXTC:
if (fr->bX)
{
snew(xout, nind);
for (i = 0; i < nind; i++)
{
copy_rvec(fr->x[ind[i]], xout[i]);
}
}
break;
default:
break;
}
switch (ftp)
{
case efTNG:
gmx_write_tng_from_trxframe(status->tng, fr, nind);
break;
case efXTC:
write_xtc(status->fio, nind, fr->step, fr->time, fr->box, xout, prec);
break;
case efTRR:
gmx_trr_write_frame(status->fio, nframes_read(status),
fr->time, fr->step, fr->box, nind, xout, vout, fout);
break;
case efGRO:
case efPDB:
case efBRK:
case efENT:
if (!fr->bAtoms)
{
gmx_fatal(FARGS, "Can not write a %s file without atom names",
ftp2ext(ftp));
}
//.........这里部分代码省略.........
开发者ID:MrTheodor,项目名称:gromacs,代码行数:101,代码来源:trxio.cpp
示例4: mk_diamond
static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab,
gmx_bool bPBC,matrix box)
{
int i,ib,j,k,l,m,nrm=0;
t_bbb *bbb;
gmx_bool *bRemove;
rvec dx;
do {
nrm = 0;
bbb = mk_bonds(a->nr,x,odist,bPBC,box);
for(i=0; (i<a->nr); i++) {
if (bbb[i].n < 2) {
for(k=0; (k<bbb[i].n); k++) {
ib = bbb[i].aa[k];
for(j=0; (j<bbb[ib].n); j++)
if (bbb[ib].aa[j] == i)
break;
if (j == bbb[ib].n)
gmx_fatal(FARGS,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
for( ; (j<bbb[ib].n-1); j++)
bbb[ib].aa[j] = bbb[ib].aa[j+1];
bbb[ib].n--;
nrm++;
}
bbb[i].n = 0;
}
}
for(i=j=0; (i<a->nr); i++) {
if (bbb[i].n >= 2) {
copy_rvec(x[i],x[j]);
j++;
}
}
fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
a->nr-j,a->nr);
a->nr = j;
sfree(bbb);
} while (nrm > 0);
/* Rename atoms */
bbb = mk_bonds(a->nr,x,odist,bPBC,box);
for(i=0; (i<a->nr); i++) {
switch (bbb[i].n) {
case 4:
a->atomname[i] = put_symtab(symtab,"C");
break;
case 3:
a->atomname[i] = put_symtab(symtab,"CH1");
break;
case 2:
a->atomname[i] = put_symtab(symtab,"CH2");
break;
default:
gmx_fatal(FARGS,"This atom (%d) has %d bonds only",i,bbb[i].n);
}
}
sfree(bbb);
}
开发者ID:daniellandau,项目名称:gromacs,代码行数:61,代码来源:mkice.c
示例5: gmx_bundle
//.........这里部分代码省略.........
xvgr_tlabel(),"(degrees)");
if (bPrintXvgrCodes())
fprintf(fkinkr,"@ subtitle \"+ = ) ( - = ( )\"\n");
fkinkl = xvgropen(opt2fn("-okl",NFILE,fnm),"Lateral kink angles",
xvgr_tlabel(),"(degrees)");
}
if (opt2bSet("-oa",NFILE,fnm)) {
init_t_atoms(&outatoms,3*n,FALSE);
outatoms.nr = 3*n;
for(i=0; i<3*n; i++) {
outatoms.atomname[i] = &anm;
outatoms.atom[i].resnr = i/3;
outatoms.resname[i/3] = &rnm;
}
fpdb = open_trx(opt2fn("-oa",NFILE,fnm),"w");
} else
fpdb = -1;
read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_NEED_X);
do {
rm_pbc(&top.idef,ePBC,fr.natoms,fr.box,fr.x,fr.x);
calc_axes(fr.x,top.atoms.atom,gnx,index,!bZ,&bun);
t = convert_time(fr.time);
fprintf(flen," %10g",t);
fprintf(fdist," %10g",t);
fprintf(fz," %10g",t);
fprintf(ftilt," %10g",t);
fprintf(ftiltr," %10g",t);
fprintf(ftiltl," %10g",t);
if (bKink) {
fprintf(fkink," %10g",t);
fprintf(fkinkr," %10g",t);
fprintf(fkinkl," %10g",t);
}
for(i=0; i<bun.n; i++) {
fprintf(flen," %6g",bun.len[i]);
fprintf(fdist," %6g",norm(bun.mid[i]));
fprintf(fz," %6g",bun.mid[i][ZZ]);
fprintf(ftilt," %6g",RAD2DEG*acos(bun.dir[i][ZZ]));
comp = bun.mid[i][XX]*bun.dir[i][XX]+bun.mid[i][YY]*bun.dir[i][YY];
fprintf(ftiltr," %6g",RAD2DEG*
asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
comp = bun.mid[i][YY]*bun.dir[i][XX]-bun.mid[i][XX]*bun.dir[i][YY];
fprintf(ftiltl," %6g",RAD2DEG*
asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
if (bKink) {
rvec_sub(bun.end[0][i],bun.end[2][i],va);
rvec_sub(bun.end[2][i],bun.end[1][i],vb);
unitv_no_table(va,va);
unitv_no_table(vb,vb);
fprintf(fkink," %6g",RAD2DEG*acos(iprod(va,vb)));
cprod(va,vb,vc);
copy_rvec(bun.mid[i],vr);
vr[ZZ] = 0;
unitv_no_table(vr,vr);
fprintf(fkinkr," %6g",RAD2DEG*asin(iprod(vc,vr)));
vl[XX] = vr[YY];
vl[YY] = -vr[XX];
vl[ZZ] = 0;
fprintf(fkinkl," %6g",RAD2DEG*asin(iprod(vc,vl)));
}
}
fprintf(flen,"\n");
fprintf(fdist,"\n");
fprintf(fz,"\n");
fprintf(ftilt,"\n");
fprintf(ftiltr,"\n");
fprintf(ftiltl,"\n");
if (bKink) {
fprintf(fkink,"\n");
fprintf(fkinkr,"\n");
fprintf(fkinkl,"\n");
}
if (fpdb >= 0)
dump_axes(fpdb,&fr,&outatoms,&bun);
} while(read_next_frame(status,&fr));
close_trx(status);
if (fpdb >= 0)
close_trx(fpdb);
fclose(flen);
fclose(fdist);
fclose(fz);
fclose(ftilt);
fclose(ftiltr);
fclose(ftiltl);
if (bKink) {
fclose(fkink);
fclose(fkinkr);
fclose(fkinkl);
}
thanx(stderr);
return 0;
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:101,代码来源:gmx_bundle.c
示例6: read_eigenvectors
void read_eigenvectors(const char *file, int *natoms, gmx_bool *bFit,
rvec **xref, gmx_bool *bDMR,
rvec **xav, gmx_bool *bDMA,
int *nvec, int **eignr,
rvec ***eigvec, real **eigval)
{
gmx_trr_header_t head;
int i, snew_size;
struct t_fileio *status;
rvec *x;
matrix box;
gmx_bool bOK;
*bDMR = FALSE;
/* read (reference (t=-1) and) average (t=0) structure */
status = gmx_trr_open(file, "r");
gmx_trr_read_frame_header(status, &head, &bOK);
*natoms = head.natoms;
snew(*xav, *natoms);
gmx_trr_read_frame_data(status, &head, box, *xav, NULL, NULL);
if ((head.t >= -1.1) && (head.t <= -0.9))
{
snew(*xref, *natoms);
for (i = 0; i < *natoms; i++)
{
copy_rvec((*xav)[i], (*xref)[i]);
}
*bDMR = (head.lambda > 0.5);
*bFit = (head.lambda > -0.5);
if (*bFit)
{
fprintf(stderr, "Read %smass weighted reference structure with %d atoms from %s\n", *bDMR ? "" : "non ", *natoms, file);
}
else
{
fprintf(stderr, "Eigenvectors in %s were determined without fitting\n", file);
sfree(*xref);
*xref = NULL;
}
gmx_trr_read_frame_header(status, &head, &bOK);
gmx_trr_read_frame_data(status, &head, box, *xav, NULL, NULL);
}
else
{
*bFit = TRUE;
*xref = NULL;
}
*bDMA = (head.lambda > 0.5);
if ((head.t <= -0.01) || (head.t >= 0.01))
{
fprintf(stderr, "WARNING: %s does not start with t=0, which should be the "
"average structure. This might not be a eigenvector file. "
"Some things might go wrong.\n",
file);
}
else
{
fprintf(stderr,
"Read %smass weighted average/minimum structure with %d atoms from %s\n",
*bDMA ? "" : "non ", *natoms, file);
}
snew(x, *natoms);
snew_size = 10;
snew(*eignr, snew_size);
snew(*eigval, snew_size);
snew(*eigvec, snew_size);
*nvec = 0;
while (gmx_trr_read_frame_header(status, &head, &bOK))
{
gmx_trr_read_frame_data(status, &head, box, x, NULL, NULL);
if (*nvec >= snew_size)
{
snew_size += 10;
srenew(*eignr, snew_size);
srenew(*eigval, snew_size);
srenew(*eigvec, snew_size);
}
i = head.step;
(*eigval)[*nvec] = head.t;
(*eignr)[*nvec] = i-1;
snew((*eigvec)[*nvec], *natoms);
for (i = 0; i < *natoms; i++)
{
copy_rvec(x[i], (*eigvec)[*nvec][i]);
}
(*nvec)++;
}
sfree(x);
gmx_trr_close(status);
fprintf(stderr, "Read %d eigenvectors (for %d atoms)\n\n", *nvec, *natoms);
}
开发者ID:MichalKononenko,项目名称:gromacs,代码行数:95,代码来源:eigio.cpp
示例7: calc_running_com
void calc_running_com(t_pull *pull) {
int i,j,n;
rvec ave;
real tm;
/* for group i, we have nhist[i] points of history. The maximum nr of points
is pull->reflag. The array comhist[i][0..nhist[i]-1] has the positions
of the center of mass over the last nhist[i] steps. x_unc[i] has the
new ones. Remove the oldest one from the list, add the new one, calculate
the average, put that in x_unc instead and return. We just move all coms
1 down in the list and add the latest one to the top.
*/
if (pull->bCyl) {
/* act on dyna groups */
for (i=0;i<pull->pull.n;i++) {
clear_rvec(ave);
for (j=0;j<(pull->reflag-1);j++) {
copy_rvec(pull->dyna.comhist[i][j+1],pull->dyna.comhist[i][j]);
rvec_add(ave,pull->dyna.comhist[i][j],ave);
}
copy_rvec(pull->dyna.x_unc[i],pull->dyna.comhist[i][j]);
rvec_add(ave,pull->dyna.x_unc[i],ave);
svmul(1.0/pull->reflag,ave,ave);
/* now ave has the running com for group i, copy it to x_unc[i] */
copy_rvec(ave,pull->dyna.x_unc[i]);
#ifdef DEBUG
if (pull->bVerbose)
for (n=0;n<pull->reflag;n++)
fprintf(stderr,"Comhist %d, grp %d: %8.3f%8.3f%8.3f\n",n,i,
pull->dyna.comhist[i][n][XX],
pull->dyna.comhist[i][n][YY],
pull->dyna.comhist[i][n][ZZ]);
#endif
}
} else {
/* act on ref group */
clear_rvec(ave);
for (j=0;j<(pull->reflag-1);j++) {
copy_rvec(pull->ref.comhist[0][j+1],pull->ref.comhist[0][j]);
rvec_add(ave,pull->ref.comhist[0][j],ave);
}
copy_rvec(pull->ref.x_unc[0],pull->ref.comhist[0][j]);
rvec_add(ave,pull->ref.x_unc[0],ave);
svmul(1.0/pull->reflag,ave,ave);
/* now ave has the running com for group i, copy it to x_unc[0] */
copy_rvec(ave,pull->ref.x_unc[0]);
#ifdef DEBUG
if (pull->bVerbose)
for (i=0;i<pull->reflag;i++)
fprintf(stderr,"Comhist %d: %8.3f%8.3f%8.3f\n",i,
pull->ref.comhist[0][i][XX],
pull->ref.comhist[0][i][YY],
pull->ref.comhist[0][i][ZZ]);
#endif
}
}
开发者ID:Chadi-akel,项目名称:cere,代码行数:62,代码来源:pullutil.c
示例8: dd_pmeredist_pos_coeffs
static void dd_pmeredist_pos_coeffs(struct gmx_pme_t *pme,
int n, gmx_bool bX, rvec *x, real *data,
pme_atomcomm_t *atc)
{
int *commnode, *buf_index;
int nnodes_comm, i, nsend, local_pos, buf_pos, node, scount, rcount;
commnode = atc->node_dest;
buf_index = atc->buf_index;
nnodes_comm = min(2*atc->maxshift, atc->nslab-1);
nsend = 0;
for (i = 0; i < nnodes_comm; i++)
{
buf_index[commnode[i]] = nsend;
nsend += atc->count[commnode[i]];
}
if (bX)
{
if (atc->count[atc->nodeid] + nsend != n)
{
gmx_fatal(FARGS, "%d particles communicated to PME rank %d are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension %c.\n"
"This usually means that your system is not well equilibrated.",
n - (atc->count[atc->nodeid] + nsend),
pme->nodeid, 'x'+atc->dimind);
}
if (nsend > pme->buf_nalloc)
{
pme->buf_nalloc = over_alloc_dd(nsend);
srenew(pme->bufv, pme->buf_nalloc);
srenew(pme->bufr, pme->buf_nalloc);
}
atc->n = atc->count[atc->nodeid];
for (i = 0; i < nnodes_comm; i++)
{
scount = atc->count[commnode[i]];
/* Communicate the count */
if (debug)
{
fprintf(debug, "dimind %d PME rank %d send to rank %d: %d\n",
atc->dimind, atc->nodeid, commnode[i], scount);
}
pme_dd_sendrecv(atc, FALSE, i,
&scount, sizeof(int),
&atc->rcount[i], sizeof(int));
atc->n += atc->rcount[i];
}
pme_realloc_atomcomm_things(atc);
}
local_pos = 0;
for (i = 0; i < n; i++)
{
node = atc->pd[i];
if (node == atc->nodeid)
{
/* Copy direct to the receive buffer */
if (bX)
{
copy_rvec(x[i], atc->x[local_pos]);
}
atc->coefficient[local_pos] = data[i];
local_pos++;
}
else
{
/* Copy to the send buffer */
if (bX)
{
copy_rvec(x[i], pme->bufv[buf_index[node]]);
}
pme->bufr[buf_index[node]] = data[i];
buf_index[node]++;
}
}
buf_pos = 0;
for (i = 0; i < nnodes_comm; i++)
{
scount = atc->count[commnode[i]];
rcount = atc->rcount[i];
if (scount > 0 || rcount > 0)
{
if (bX)
{
/* Communicate the coordinates */
pme_dd_sendrecv(atc, FALSE, i,
pme->bufv[buf_pos], scount*sizeof(rvec),
atc->x[local_pos], rcount*sizeof(rvec));
}
/* Communicate the coefficients */
pme_dd_sendrecv(atc, FALSE, i,
pme->bufr+buf_pos, scount*sizeof(real),
atc->coefficient+local_pos, rcount*sizeof(real));
buf_pos += scount;
local_pos += atc->rcount[i];
//.........这里部分代码省略.........
开发者ID:FoldingAtHome,项目名称:gromacs,代码行数:101,代码来源:pme-redistribute.c
示例9: dd_pmeredist_f
void dd_pmeredist_f(struct gmx_pme_t *pme, pme_atomcomm_t *atc,
int n, rvec *f,
gmx_bool bAddF)
{
int *commnode, *buf_index;
int nnodes_comm, local_pos, buf_pos, i, scount, rcount, node;
commnode = atc->node_dest;
buf_index = atc->buf_index;
nnodes_comm = min(2*atc->maxshift, atc->nslab-1);
local_pos = atc->count[atc->nodeid];
buf_pos = 0;
for (i = 0; i < nnodes_comm; i++)
{
scount = atc->rcount[i];
rcount = atc->count[commnode[i]];
if (scount > 0 || rcount > 0)
{
/* Communicate the forces */
pme_dd_sendrecv(atc, TRUE, i,
atc->f[local_pos], scount*sizeof(rvec),
pme->bufv[buf_pos], rcount*sizeof(rvec));
local_pos += scount;
}
buf_index[commnode[i]] = buf_pos;
buf_pos += rcount;
}
local_pos = 0;
if (bAddF)
{
for (i = 0; i < n; i++)
{
node = atc->pd[i];
if (node == atc->nodeid)
{
/* Add from the local force array */
rvec_inc(f[i], atc->f[local_pos]);
local_pos++;
}
else
{
/* Add from the receive buffer */
rvec_inc(f[i], pme->bufv[buf_index[node]]);
buf_index[node]++;
}
}
}
else
{
for (i = 0; i < n; i++)
{
node = atc->pd[i];
if (node == atc->nodeid)
{
/* Copy from the local force array */
copy_rvec(atc->f[local_pos], f[i]);
local_pos++;
}
else
{
/* Copy from the receive buffer */
copy_rvec(pme->bufv[buf_index[node]], f[i]);
buf_index[node]++;
}
}
}
}
开发者ID:FoldingAtHome,项目名称:gromacs,代码行数:70,代码来源:pme-redistribute.c
示例10: gmx_helixorient
//.........这里部分代码省略.........
fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
"Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
oenv);
fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
"Incremental local helix rotation", "Time(ps)",
"Rotation (degrees)", oenv);
}
else
{
fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
"Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
"Cumulative local helix rotation", "Time(ps)",
"Rotation (degrees)", oenv);
}
clear_rvecs(3, unitaxes);
unitaxes[0][0] = 1;
unitaxes[1][1] = 1;
unitaxes[2][2] = 1;
gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
do
{
/* initialisation for correct distance calculations */
set_pbc(&pbc, ePBC, box);
/* make molecules whole again */
gmx_rmpbc(gpbc, natoms, box, x);
/* copy coords to our smaller arrays */
for (i = 0; i < iCA; i++)
{
copy_rvec(x[ind_CA[i]], x_CA[i]);
if (bSC)
{
copy_rvec(x[ind_SC[i]], x_SC[i]);
}
}
for (i = 0; i < iCA-3; i++)
{
rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
rvec_sub(r12[i], r23[i], diff13[i]);
rvec_sub(r23[i], r34[i], diff24[i]);
/* calculate helix axis */
cprod(diff13[i], diff24[i], helixaxis[i]);
svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
tmp = cos_angle(diff13[i], diff24[i]);
twist[i] = 180.0/M_PI * std::acos( tmp );
radius[i] = std::sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
rise[i] = std::abs(iprod(r23[i], helixaxis[i]));
svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
}
residueradius[0] = residuetwist[0] = residuerise[0] = 0;
residueradius[1] = radius[0];
residuetwist[1] = twist[0];
开发者ID:MelroLeandro,项目名称:gromacs,代码行数:67,代码来源:gmx_helixorient.cpp
示例11: gmx_traj
//.........这里部分代码省略.........
if (bCV)
{
snew(sumv, fr.natoms);
}
if (bCF)
{
snew(sumf, fr.natoms);
}
nr_xfr = 0;
nr_vfr = 0;
nr_ffr = 0;
if (bCom && bPBC)
{
gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
}
do
{
time = output_env_conv_time(oenv, fr.time);
if (fr.bX && bNoJump && fr.bBox)
{
if (xp)
{
remove_jump(fr.box, fr.natoms, xp, fr.x);
}
else
{
snew(xp, fr.natoms);
}
for (i = 0; i < fr.natoms; i++)
{
copy_rvec(fr.x[i], xp[i]);
}
}
if (fr.bX && bCom && bPBC)
{
gmx_rmpbc_trxfr(gpbc, &fr);
}
if (bVD && fr.bV)
{
update_histo(isize[0], index[0], fr.v, &nvhisto, &vhisto, binwidth);
}
if (bOX && fr.bX)
{
print_data(outx, time, fr.x, mass, bCom, ngroups, isize, index, bDim, sffmt);
}
if (bOXT && fr.bX)
{
frout = fr;
if (!frout.bAtoms)
{
frout.atoms = &top.atoms;
frout.bAtoms = TRUE;
}
write_trx_x(status_out, &frout, mass, bCom, ngroups, isize, index);
}
if (bOV && fr.bV)
{
print_data(outv, time, fr.v, mass, bCom, ngroups, isize, index, bDim, sffmt);
}
if (bOF && fr.bF)
开发者ID:MelroLeandro,项目名称:gromacs,代码行数:67,代码来源:gmx_traj.cpp
示例12: connelly_plot
void connelly_plot(const char *fn,int ndots,real dots[],rvec x[],t_atoms *atoms,
t_symtab *symtab,int ePBC,matrix box,gmx_bool bSave)
{
static const char *atomnm="DOT";
static const char *resnm ="DOT";
static const char *title ="Connely Dot Surface Generated by g_sas";
int i,i0,r0,ii0,k;
rvec *xnew;
t_atoms aaa;
if (bSave) {
i0 = atoms->nr;
r0 = atoms->nres;
srenew(atoms->atom,atoms->nr+ndots);
srenew(atoms->atomname,atoms->nr+ndots);
srenew(atoms->resinfo,r0+1);
atoms->atom[i0].resind = r0;
t_atoms_set_resinfo(atoms,i0,symtab,resnm,r0+1,' ',0,' ');
srenew(atoms->pdbinfo,atoms->nr+ndots);
snew(xnew,atoms->nr+ndots);
for(i=0; (i<atoms->nr); i++)
copy_rvec(x[i],xnew[i]);
for(i=k=0; (i<ndots); i++) {
ii0 = i0+i;
atoms->atomname[ii0] = put_symtab(symtab,atomnm);
atoms->pdbinfo[ii0].type = epdbATOM;
atoms->pdbinfo[ii0].atomnr= ii0;
atoms->atom[ii0].resind = r0;
xnew[ii0][XX] = dots[k++];
xnew[ii0][YY] = dots[k++];
xnew[ii0][ZZ] = dots[k++];
atoms->pdbinfo[ii0].bfac = 0.0;
atoms->pdbinfo[ii0].occup = 0.0;
}
atoms->nr = i0+ndots;
atoms->nres = r0+1;
write_sto_conf(fn,title,atoms,xnew,NULL,ePBC,box);
atoms->nres = r0;
atoms->nr = i0;
}
else {
init_t_atoms(&aaa,ndots,TRUE);
aaa.atom[0].resind = 0;
t_atoms_set_resinfo(&aaa,0,symtab,resnm,1,' ',0,' ');
snew(xnew,ndots);
for(i=k=0; (i<ndots); i++) {
ii0 = i;
aaa.atomname[ii0] = put_symtab(symtab,atomnm);
aaa.pdbinfo[ii0].type = epdbATOM;
aaa.pdbinfo[ii0].atomnr= ii0;
aaa.atom[ii0].resind = 0;
xnew[ii0][XX] = dots[k++];
xnew[ii0][YY] = dots[k++];
xnew[ii0][ZZ] = dots[k++];
aaa.pdbinfo[ii0].bfac = 0.0;
aaa.pdbinfo[ii0].occup = 0.0;
}
aaa.nr = ndots;
write_sto_conf(fn,title,&aaa,xnew,NULL,ePBC,box);
do_conect(fn,ndots,xnew);
free_t_atoms(&aaa,FALSE);
}
sfree(xnew);
}
开发者ID:BradleyDickson,项目名称:ABPenabledGROMACS,代码行数:65,代码来源:gmx_sas.c
示例13: sgangle_plot_single
void sgangle_plot_single(const char *fn,const char *afile,const char *dfile,
const char *d1file, const char *d2file,
atom_id index1[], int gnx1, char *grpn1,
atom_id index2[], int gnx2, char *grpn2,
t_topology *top,int ePBC, const output_env_t oenv)
{
FILE
*sg_angle, /* xvgr file with angles */
*sg_distance = NULL, /* xvgr file with distances */
*sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
*sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
real
t, /* time */
angle, /* cosine of angle between two groups */
distance, /* distance between two groups. */
distance1, /* distance between plane and one of two atoms */
distance2; /* same for second of two atoms */
t_trxstatus *status;
int natoms,teller=0;
int i;
rvec *x0; /* coordinates, and coordinates corrected for pb */
rvec *xzero;
matrix box;
char buf[256]; /* for xvgr title */
gmx_rmpbc_t gpbc=NULL;
if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
if (dfile) {
sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
}
if (d1file) {
sprintf(buf,"Distance between plane and first atom of vector");
sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
}
if (d2file) {
sprintf(buf,"Distance between plane and second atom of vector");
sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
}
snew(xzero,natoms);
gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
do {
teller++;
gmx_rmpbc(gpbc,natoms,box,x0);
if (teller==1) {
for(i=0;i<natoms;i++)
copy_rvec(x0[i],xzero[i]);
}
calc_angle_single(ePBC,box,
xzero,x0,index1,index2,gnx1,gnx2,&angle,
&distance,&distance1,&distance2);
fprintf(sg_angle,"%12g %12g %12g\n",t,angle,acos(angle)*180.0/M_PI);
if (dfile)
fprintf(sg_distance,"%12g %12g\n",t,distance);
if (d1file)
fprintf(sg_distance1,"%12g %12g\n",t,distance1);
if (d2file)
fprintf(sg_distance2,"%12g %12g\n",t,distance1);
} while (read_next_x(oenv,status,&t,natoms,x0,box));
gmx_rmpbc_done(gpbc);
fprintf(stderr,"\n");
close_trj(status);
ffclose(sg_angle);
if (dfile)
ffclose(sg_distance);
if (d1file)
ffclose(sg_distance1);
if (d2file)
ffclose(sg_distance2);
sfree(x0);
}
开发者ID:martinhoefling,项目名称:gromacs,代码行数:88,代码来源:gmx_sgangle.c
示例14: pull_calc_coms
/* calculates center of mass of selection index from all coordinates x */
void pull_calc_coms(t_commrec *cr,
struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, double t,
rvec x[], rvec *xp)
{
int g;
real twopi_box = 0;
pull_comm_t *comm;
comm = &pull->comm;
if (comm->rbuf == NULL)
{
snew(comm->rbuf, pull->ngroup);
}
if (comm->dbuf == NULL)
{
snew(comm->dbuf, 3*pull->ngroup);
}
if (pull->bRefAt && pull->bSetPBCatoms)
{
pull_set_pbcatoms(cr, pull, x, comm->rbuf);
if (cr != NULL && DOMAINDECOMP(cr))
{
/* We can keep these PBC reference coordinates fixed for nstlist
* steps, since atoms won't jump over PBC.
* This avoids a global reduction at the next nstlist-1 steps.
* Note that the exact values of the pbc reference coordinates
* are irrelevant, as long all atoms in the group are within
* half a box distance of the reference coordinate.
*/
pull->bSetPBCatoms = FALSE;
}
}
if (pull->cosdim >= 0)
{
int m;
assert(pull->npbcdim <= DIM);
for (m = pull->cosdim+1; m < pull->npbcdim; m++)
{
if (pbc->box[m][pull->cosdim] != 0)
{
gmx_fatal(FARGS, "Can not do cosine weighting for trilinic dimensions");
}
}
twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
}
for (g = 0; g < pull->ngroup; g++)
{
pull_group_work_t *pgrp;
pgrp = &pull->group[g];
if (pgrp->bCalcCOM)
{
if (pgrp->epgrppbc != epgrppbcCOS)
{
dvec com, comp;
double wmass, wwmass;
rvec x_pbc = { 0, 0, 0 };
int i;
clear_dvec(com);
clear_dvec(comp);
wmass = 0;
wwmass = 0;
if (pgrp->epgrppbc == epgrppbcREFAT)
{
/* Set the pbc atom */
copy_rvec(comm->rbuf[g], x_pbc);
}
for (i = 0; i < pgrp->nat_loc; i++)
{
int ii, m;
real mass, wm;
ii = pgrp->ind_loc[i];
mass = md->massT[ii];
if (pgrp->weight_loc == NULL)
{
wm = mass;
wmass += wm;
}
else
{
real w;
w = pgrp->weight_loc[i];
wm = w*mass;
wmass += wm;
wwmass += wm*w;
}
//.........这里部分代码省略.........
开发者ID:MelroLeandro,项目名称:gromacs,代码行数:101,代码来源:pullutil.cpp
示例15: calc_h_pos
//.........这里部分代码省略.........
break;
case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
for (d = 0; (d < DIM); d++)
{
xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
xH2[d] = ( xAI[d]
- distH*sin(alfaH)*0.5*sb[d]
+ distH*sin(alfaH)*s6*sa[d]
- distH*cos(alfaH)*sij[d] );
if (xH3[XX] != NOTSET && xH3[YY] != NOTSET && xH3[ZZ] != NOTSET)
{
xH3[d] = ( xAI[d]
- distH*sin(alfaH)*0.5*sb[d]
- distH*sin(alfaH)*s6*sa[d]
- distH*cos(alfaH)*sij[d] );
}
}
break;
case 5: /* one tetrahedral hydrogen, e.g. C3CH */
{
real center;
rvec dxc;
for (d = 0; (d < DIM); d++)
{
center = (xAJ[d]+xAK[d]+xAL[d])/3.0;
dxc[d] = xAI[d]-center;
}
center = norm(dxc);
for (d = 0; (d < DIM); d++)
{
xH1[d] = xAI[d]+dxc[d]*distH/center;
}
break;
}
case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
{
rvec rBB, rCC1, rCC2, rNN;
real bb, nn;
for (d = 0; (d < DIM); d++)
{
rBB[d] = xAI[d]-0.5*(xAJ[d]+xAK[d]);
}
bb = norm(rBB);
rvec_sub(xAI, xAJ, rCC1);
rvec_sub(xAI, xAK, rCC2);
cprod(rCC1, rCC2, rNN);
nn = norm(rNN);
for (d = 0; (d < DIM); d++)
{
xH1[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
sin(alfaH/2.0)*rNN[d]/nn);
xH2[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
sin(alfaH/2.0)*rNN[d]/nn);
}
break;
}
case 7: /* two water hydrogens */
gen_waterhydrogen(2, xa, xh, l);
break;
case 10: /* three water hydrogens */
gen_waterhydrogen(3, xa, xh, l);
break;
case 11: /* four water hydrogens */
gen_waterhydrogen(4, xa, xh, l);
break;
case 8: /* two carboxyl oxygens, -COO- */
for (d = 0; (d < DIM); d++)
{
xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
}
break;
case 9: /* carboxyl oxygens and hydrogen, -COOH */
{
rvec xa2[4]; /* i,j,k,l */
/* first add two oxygens */
for (d = 0; (d < DIM); d++)
{
xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
}
/* now use rule 2 to add hydrogen to 2nd oxygen */
copy_rvec(xH2, xa2[0]); /* new i = n' */
copy_rvec(xAI, xa2[1]); /* new j = i */
copy_rvec(xAJ, xa2[2]); /* new k = j */
copy_rvec(xAK, xa2[3]); /* new l = k, not used */
calc_h_pos(2, xa2, (xh+2), l);
break;
}
default:
gmx_fatal(FARGS, "Invalid argument (%d) for nht in routine genh\n", nht);
} /* end switch */
}
开发者ID:yhalcyon,项目名称:Gromacs,代码行数:101,代码来源:calch.c
示例16: dielectric
//.........这里部分代码省略.........
GMX_RELEASE_ASSERT(time != NULL, "Memory not allocated correctly - time array is NULL");
if (nfr == 0)
{
t0 = fr.time;
}
time[nfr] = fr.time-t0;
if (time[nfr] <= bfit)
{
bi = nfr;
}
if (time[nfr] <= efit)
{
ei = nfr;
}
if (bNoJump)
{
if (xp)
{
remove_jump(fr.box, fr.natoms, xp, fr.x);
}
else
{
snew(xp, fr.natoms);
}
for (i = 0; i < fr.natoms; i++)
{
copy_rvec(fr.x[i], xp[i]);
}
}
gmx_rmpbc_trxfr(gpbc, &fr);
calc_mj(top, ePBC, fr.box, bNoJump, nmols, indexm, fr.x, mtrans[nfr], mass2, qmol);
for (i = 0; i < isize; i++)
{
j = index0[i];
svmul(top.atoms.atom[j].q, fr.x[j], fr.x[j]);
rvec_inc(mu[nfr], fr.x[j]);
}
/*if(mod(nfr,nshift)==0){*/
if (nfr%nshift == 0)
{
for (j = nfr; j >= 0; j--)
{
rvec_sub(mtrans[nfr], mtrans[j], tmp);
dsp2[nfr-j] += norm2(tmp);
xshfr[nfr-j] += 1.0;
}
}
if (fr.bV)
{
if (nvfr >= valloc)
{
valloc += 100;
srenew(vfr, valloc);
开发者ID:pjohansson,项目名称:gromacs,代码行数:67,代码来源:gmx_current.cpp
|
请发表评论