本文整理汇总了C++中set_pbc函数的典型用法代码示例。如果您正苦于以下问题：C++ set_pbc函数的具体用法？C++ set_pbc怎么用？C++ set_pbc使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了set_pbc函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: set_pbc
static t_bbb *mk_bonds(int natoms,rvec x[],real odist,
gmx_bool bPBC,matrix box)
{
real od2 = odist*odist+1e-5;
t_pbc pbc;
t_bbb *bbb;
int i,j;
rvec dx;
if (bPBC)
set_pbc(&pbc,box);
snew(bbb,natoms);
for(i=0; (i<natoms); i++) {
for(j=i+1; (j<natoms); j++) {
if (bPBC)
pbc_dx(&pbc,x[i],x[j],dx);
else
rvec_sub(x[i],x[j],dx);
if (iprod(dx,dx) <= od2) {
bbb[i].aa[bbb[i].n++] = j;
bbb[j].aa[bbb[j].n++] = i;
}
}
}
if (debug)
#define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
for(i=0; (i<natoms); i++)
fprintf(debug,"bonds from %3d: %d %d %d %d\n",
i,PRB(1),PRB(2),PRB(3),PRB(4));
#undef PRB
return bbb;
}
开发者ID:daniellandau，项目名称:gromacs，代码行数:32，代码来源:mkice.c
示例2: allPairsDistOk
/* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
* leaks memory (May 2012). The function could be deleted as soon as the momory leaks
* in addconf.c are fixed.
* However, when inserting a small molecule in a system containing not too many atoms,
* allPairsDistOk is probably even faster than addconf.c
*/
static gmx_bool
allPairsDistOk(t_atoms *atoms, rvec *x, real *r,
int ePBC, matrix box,
t_atoms *atoms_insrt, rvec *x_n, real *r_insrt)
{
int i, j;
rvec dx;
real n2, r2;
t_pbc pbc;
set_pbc(&pbc, ePBC, box);
for (i = 0; i < atoms->nr; i++)
{
for (j = 0; j < atoms_insrt->nr; j++)
{
pbc_dx(&pbc, x[i], x_n[j], dx);
n2 = norm2(dx);
r2 = sqr(r[i]+r_insrt[j]);
if (n2 < r2)
{
return FALSE;
}
}
}
return TRUE;
}
开发者ID:alwanderer，项目名称:gromacs，代码行数:32，代码来源:gmx_genbox.cpp
示例3: set_pull_init
void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,
const output_env_t oenv)
{
pull_params_t *pull;
struct pull_t *pull_work;
t_mdatoms *md;
t_pbc pbc;
int c;
double t_start;
pull = ir->pull;
pull_work = init_pull(NULL, pull, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
md = init_mdatoms(NULL, mtop, ir->efep);
atoms2md(mtop, ir, 0, NULL, mtop->natoms, md);
if (ir->efep)
{
update_mdatoms(md, lambda);
}
set_pbc(&pbc, ir->ePBC, box);
t_start = ir->init_t + ir->init_step*ir->delta_t;
pull_calc_coms(NULL, pull_work, md, &pbc, t_start, x, NULL);
fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
for (c = 0; c < pull->ncoord; c++)
{
t_pull_coord *pcrd;
t_pull_group *pgrp0, *pgrp1;
double value;
real init = 0;
pcrd = &pull->coord[c];
pgrp0 = &pull->group[pcrd->group[0]];
pgrp1 = &pull->group[pcrd->group[1]];
fprintf(stderr, "%8d %8d %8d\n",
pcrd->group[0], pgrp0->nat, pgrp0->pbcatom+1);
fprintf(stderr, "%8d %8d %8d ",
pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1);
if (pcrd->bStart)
{
init = pcrd->init;
pcrd->init = 0;
}
get_pull_coord_value(pull_work, c, &pbc, &value);
fprintf(stderr, " %10.3f nm", value);
if (pcrd->bStart)
{
pcrd->init = value + init;
}
fprintf(stderr, " %10.3f nm\n", pcrd->init);
}
finish_pull(pull_work);
}
开发者ID:JehandadKhan，项目名称:gromacs，代码行数:60，代码来源:readpull.c
示例4: calc_dist_tot
static void calc_dist_tot(int nind, atom_id index[], rvec x[],
int ePBC, matrix box,
real **d, real **dtot, real **dtot2,
gmx_bool bNMR, real **dtot1_3, real **dtot1_6)
{
int i, j;
real *xi;
real temp, temp2, temp1_3;
rvec dx;
t_pbc pbc;
set_pbc(&pbc, ePBC, box);
for (i = 0; (i < nind-1); i++)
{
xi = x[index[i]];
for (j = i+1; (j < nind); j++)
{
pbc_dx(&pbc, xi, x[index[j]], dx);
temp2 = norm2(dx);
temp = std::sqrt(temp2);
d[i][j] = temp;
dtot[i][j] += temp;
dtot2[i][j] += temp2;
if (bNMR)
{
temp1_3 = 1.0/(temp*temp2);
dtot1_3[i][j] += temp1_3;
dtot1_6[i][j] += temp1_3*temp1_3;
}
}
}
}
开发者ID:mpharrigan，项目名称:gromacs，代码行数:32，代码来源:gmx_rmsdist.cpp
示例5: mk_bonds
void mk_bonds(int nnm, t_nm2type nmt[],
t_atoms *atoms, rvec x[], t_params *bond, int nbond[],
gmx_bool bPBC, matrix box)
{
t_param b;
int i, j;
t_pbc pbc;
rvec dx;
real dx2;
for (i = 0; (i < MAXATOMLIST); i++)
{
b.a[i] = -1;
}
for (i = 0; (i < MAXFORCEPARAM); i++)
{
b.c[i] = 0.0;
}
if (bPBC)
{
set_pbc(&pbc, -1, box);
}
for (i = 0; (i < atoms->nr); i++)
{
if ((i % 10) == 0)
{
fprintf(stderr, "\ratom %d", i);
}
for (j = i+1; (j < atoms->nr); j++)
{
if (bPBC)
{
pbc_dx(&pbc, x[i], x[j], dx);
}
else
{
rvec_sub(x[i], x[j], dx);
}
dx2 = iprod(dx, dx);
if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
sqrt(dx2)))
{
b.AI = i;
b.AJ = j;
b.C0 = sqrt(dx2);
add_param_to_list (bond, &b);
nbond[i]++;
nbond[j]++;
if (debug)
{
fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
*atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
}
}
}
}
fprintf(stderr, "\ratom %d\n", i);
}
开发者ID:exianshine，项目名称:gromacs，代码行数:60，代码来源:g_x2top.c
示例6: chk_bonds
static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
{
int ftype, i, k, ai, aj, type;
real b0, blen, deviation, devtot;
t_pbc pbc;
rvec dx;
devtot = 0;
set_pbc(&pbc, ePBC, box);
for (ftype = 0; (ftype < F_NRE); ftype++)
{
if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
{
for (k = 0; (k < idef->il[ftype].nr); )
{
type = idef->il[ftype].iatoms[k++];
ai = idef->il[ftype].iatoms[k++];
aj = idef->il[ftype].iatoms[k++];
b0 = 0;
switch (ftype)
{
case F_BONDS:
b0 = idef->iparams[type].harmonic.rA;
break;
case F_G96BONDS:
b0 = sqrt(idef->iparams[type].harmonic.rA);
break;
case F_MORSE:
b0 = idef->iparams[type].morse.b0A;
break;
case F_CUBICBONDS:
b0 = idef->iparams[type].cubic.b0;
break;
case F_CONSTR:
b0 = idef->iparams[type].constr.dA;
break;
default:
break;
}
if (b0 != 0)
{
pbc_dx(&pbc, x[ai], x[aj], dx);
blen = norm(dx);
deviation = sqr(blen-b0);
if (sqrt(deviation/sqr(b0) > tol))
{
fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
}
}
}
}
}
}
开发者ID:dkarkoulis，项目名称:gromacs，代码行数:53，代码来源:check.c
示例7: calc_mj
static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
rvec fr[], rvec mj, real mass2[], real qmol[])
{
int i, j, k, l;
rvec tmp;
rvec center;
rvec mt1, mt2;
t_pbc pbc;
calc_box_center(ecenterRECT, box, center);
if (!bNoJump)
{
set_pbc(&pbc, ePBC, box);
}
clear_rvec(tmp);
for (i = 0; i < isize; i++)
{
clear_rvec(mt1);
clear_rvec(mt2);
k = top.mols.index[index0[i]];
l = top.mols.index[index0[i+1]];
for (j = k; j < l; j++)
{
svmul(mass2[j], fr[j], tmp);
rvec_inc(mt1, tmp);
}
if (bNoJump)
{
svmul(qmol[k], mt1, mt1);
}
else
{
pbc_dx(&pbc, mt1, center, mt2);
svmul(qmol[k], mt2, mt1);
}
rvec_inc(mj, mt1);
}
}
开发者ID:pjohansson，项目名称:gromacs，代码行数:48，代码来源:gmx_current.cpp
示例8: setStateDependentAwhParams
void setStateDependentAwhParams(AwhParams *awhParams,
const pull_params_t *pull_params, pull_t *pull_work,
const matrix box, int ePBC,
const t_grpopts *inputrecGroupOptions, warninp_t wi)
{
/* The temperature is not really state depenendent but is not known
* when read_awhParams is called (in get ir).
* It is known first after do_index has been called in grompp.cpp.
*/
if (inputrecGroupOptions->ref_t == NULL ||
inputrecGroupOptions->ref_t[0] <= 0)
{
gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
}
for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
{
if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
{
gmx_fatal(FARGS, "AWH biasing is currently only supported for identical temperatures for all temperature coupling groups");
}
}
t_pbc pbc;
set_pbc(&pbc, ePBC, box);
for (int k = 0; k < awhParams->numBias; k++)
{
AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
for (int d = 0; d < awhBiasParams->ndim; d++)
{
AwhDimParams *dimParams = &awhBiasParams->dimParams[d];
/* The periodiciy of the AWH grid in certain cases depends on the simulation box */
dimParams->period = get_pull_coord_period(pull_params, dimParams->coordIndex, box);
/* The initial coordinate value, converted to external user units. */
dimParams->coordValueInit =
get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
t_pull_coord *pullCoord = &pull_params->coord[dimParams->coordIndex];
dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(pullCoord);
}
}
checkInputConsistencyInterval(awhParams, wi);
/* Register AWH as external potential with pull to check consistency. */
Awh::registerAwhWithPull(*awhParams, pull_work);
}
开发者ID:HITS-MBM，项目名称:gromacs-fda，代码行数:48，代码来源:read-params.cpp
示例9: 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
示例10: calc_mat
static void calc_mat(int nres, int natoms, int rndx[],
rvec x[], atom_id *index,
real trunc, real **mdmat, int **nmat, int ePBC, matrix box)
{
int i, j, resi, resj;
real trunc2, r, r2;
t_pbc pbc;
rvec ddx;
set_pbc(&pbc, ePBC, box);
trunc2 = sqr(trunc);
for (resi = 0; (resi < nres); resi++)
{
for (resj = 0; (resj < nres); resj++)
{
mdmat[resi][resj] = FARAWAY;
}
}
for (i = 0; (i < natoms); i++)
{
resi = rndx[i];
for (j = i+1; (j < natoms); j++)
{
resj = rndx[j];
pbc_dx(&pbc, x[index[i]], x[index[j]], ddx);
r2 = norm2(ddx);
if (r2 < trunc2)
{
nmat[resi][j]++;
nmat[resj][i]++;
}
mdmat[resi][resj] = min(r2, mdmat[resi][resj]);
}
}
for (resi = 0; (resi < nres); resi++)
{
mdmat[resi][resi] = 0;
for (resj = resi+1; (resj < nres); resj++)
{
r = sqrt(mdmat[resi][resj]);
mdmat[resi][resj] = r;
mdmat[resj][resi] = r;
}
}
}
开发者ID:yhalcyon，项目名称:Gromacs，代码行数:46，代码来源:gmx_mdmat.c
示例11: put_atoms_in_compact_unitcell
const char *
put_atoms_in_compact_unitcell(int ePBC,int ecenter,matrix box,
int natoms,rvec x[])
{
t_pbc pbc;
rvec box_center,dx;
int i;
set_pbc(&pbc,ePBC,box);
calc_box_center(ecenter,box,box_center);
for(i=0; i<natoms; i++) {
pbc_dx(&pbc,x[i],box_center,dx);
rvec_add(box_center,dx,x[i]);
}
return pbc.bLimitDistance ?
"WARNING: Could not put all atoms in the compact unitcell\n"
: NULL;
}
开发者ID:alexholehouse，项目名称:gromacs，代码行数:19，代码来源:pbc.c
示例12: calc_dist
static void calc_dist(int nind, int index[], const rvec x[], int ePBC, matrix box,
real **d)
{
int i, j;
rvec dx;
real temp2;
t_pbc pbc;
set_pbc(&pbc, ePBC, box);
for (i = 0; (i < nind-1); i++)
{
const real *xi = x[index[i]];
for (j = i+1; (j < nind); j++)
{
pbc_dx(&pbc, xi, x[index[j]], dx);
temp2 = norm2(dx);
d[i][j] = std::sqrt(temp2);
}
}
}
开发者ID:MichalKononenko，项目名称:gromacs，代码行数:20，代码来源:gmx_rmsdist.cpp
示例13: partition_absorbing_box
void
partition_absorbing_box( grid_t * g,
double gx0, double gy0, double gz0,
double gx1, double gy1, double gz1,
int gnx, int gny, int gnz,
int gpx, int gpy, int gpz,
int pbc ) {
int px, py, pz;
partition_periodic_box( g,
gx0, gy0, gz0,
gx1, gy1, gz1,
gnx, gny, gnz,
gpx, gpy, gpz );
// Override periodic boundary conditions
RANK_TO_INDEX( world_rank, px,py,pz );
if( px==0 && gnx>1 ) {
set_fbc(g,BOUNDARY(-1,0,0),absorb_fields);
set_pbc(g,BOUNDARY(-1,0,0),pbc);
}
if( px==gpx-1 && gnx>1 ) {
set_fbc(g,BOUNDARY( 1,0,0),absorb_fields);
set_pbc(g,BOUNDARY( 1,0,0),pbc);
}
if( py==0 && gny>1 ) {
set_fbc(g,BOUNDARY(0,-1,0),absorb_fields);
set_pbc(g,BOUNDARY(0,-1,0),pbc);
}
if( py==gpy-1 && gny>1 ) {
set_fbc(g,BOUNDARY(0, 1,0),absorb_fields);
set_pbc(g,BOUNDARY(0, 1,0),pbc);
}
if( pz==0 && gnz>1 ) {
set_fbc(g,BOUNDARY(0,0,-1),absorb_fields);
set_pbc(g,BOUNDARY(0,0,-1),pbc);
}
if( pz==gpz-1 && gnz>1 ) {
set_fbc(g,BOUNDARY(0,0, 1),absorb_fields);
set_pbc(g,BOUNDARY(0,0, 1),pbc);
}
}
开发者ID:Tomyao，项目名称:vpic，代码行数:49，代码来源:partition.c
示例14: put_atoms_in_compact_unitcell
void put_atoms_in_compact_unitcell(int ePBC, int ecenter, const matrix box,
int natoms, rvec x[])
{
t_pbc pbc;
rvec box_center, dx;
int i;
set_pbc(&pbc, ePBC, box);
if (pbc.ePBCDX == epbcdxUNSUPPORTED)
{
gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
}
calc_box_center(ecenter, box, box_center);
for (i = 0; i < natoms; i++)
{
pbc_dx(&pbc, x[i], box_center, dx);
rvec_add(box_center, dx, x[i]);
}
}
开发者ID:wangxubo0201，项目名称:gromacs，代码行数:21，代码来源:pbc.cpp
示例15: gmx_dist
//.........这里部分代码省略.........
if (index[g][i]>max)
max=index[g][i];
if (index[g][i] >= top->atoms.nr)
gmx_fatal(FARGS,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);
mass[g]+=top->atoms.atom[index[g][i]].m;
}
}
natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
t0 = t;
if (max>=natoms)
gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);
if (!bCutoff) {
/* open output file */
fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),
"Distance","Time (ps)","Distance (nm)",oenv);
xvgr_legend(fp,4,leg,oenv);
} else {
ngrps = 1;
if (bLifeTime)
snew(contact_time,isize[1]);
}
if (ePBC != epbcNONE)
snew(pbc,1);
else
pbc = NULL;
gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
do {
/* initialisation for correct distance calculations */
if (pbc) {
set_pbc(pbc,ePBC,box);
/* make molecules whole again */
gmx_rmpbc(gpbc,natoms,box,x);
}
/* calculate center of masses */
for(g=0;(g<ngrps);g++) {
if (isize[g] == 1) {
copy_rvec(x[index[g][0]],com[g]);
} else {
for(d=0;(d<DIM);d++) {
com[g][d]=0;
for(i=0;(i<isize[g]);i++) {
com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
}
com[g][d] /= mass[g];
}
}
}
if (!bCutoff) {
/* write to output */
fprintf(fp,"%12.7f ",t);
for(g=0;(g<ngrps/2);g++) {
if (pbc)
pbc_dx(pbc,com[2*g],com[2*g+1],dx);
else
rvec_sub(com[2*g],com[2*g+1],dx);
fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
norm(dx),dx[XX],dx[YY],dx[ZZ]);
}
fprintf(fp,"\n");
} else {
开发者ID:BradleyDickson，项目名称:ABPenabledGROMACS，代码行数:67，代码来源:gmx_dist.c
示例16: gmx_densmap
//.........这里部分代码省略.........
for(i=0; i<n1; i++)
snew(grid[i],n2);
box1 = 0;
box2 = 0;
nfr = 0;
do {
if (!bRadial) {
box1 += box[c1][c1];
box2 += box[c2][c2];
invcellvol = n1*n2;
if (nmpower == -3)
invcellvol /= det(box);
else if (nmpower == -2)
invcellvol /= box[c1][c1]*box[c2][c2];
for(i=0; i<nindex; i++) {
j = index[i];
if ((!bXmin || x[j][cav] >= xmin) &&
(!bXmax || x[j][cav] <= xmax)) {
m1 = x[j][c1]/box[c1][c1];
if (m1 >= 1)
m1 -= 1;
if (m1 < 0)
m1 += 1;
m2 = x[j][c2]/box[c2][c2];
if (m2 >= 1)
m2 -= 1;
if (m2 < 0)
m2 += 1;
grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
}
}
} else {
set_pbc(&pbc,ePBC,box);
for(i=0; i<2; i++) {
if (gnx[i] == 1) {
/* One atom, just copy the coordinates */
copy_rvec(x[ind[i][0]],xcom[i]);
} else {
/* Calculate the center of mass */
clear_rvec(xcom[i]);
mtot = 0;
for(j=0; j<gnx[i]; j++) {
k = ind[i][j];
m = top.atoms.atom[k].m;
for(l=0; l<DIM; l++)
xcom[i][l] += m*x[k][l];
mtot += m;
}
svmul(1/mtot,xcom[i],xcom[i]);
}
}
pbc_dx(&pbc,xcom[1],xcom[0],direction);
for(i=0; i<DIM; i++)
center[i] = xcom[0][i] + 0.5*direction[i];
unitv(direction,direction);
for(i=0; i<nindex; i++) {
j = index[i];
pbc_dx(&pbc,x[j],center,dx);
axial = iprod(dx,direction);
r = sqrt(norm2(dx) - axial*axial);
if (axial>=-amax && axial<amax && r<rmax) {
if (bMirror)
r += rmax;
grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
}
开发者ID:BradleyDickson，项目名称:ABPenabledGROMACS，代码行数:67，代码来源:gmx_densmap.c
示例17: read_ang_dih
void read_ang_dih(const char *trj_fn,
gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
int maxangstat, int angstat[],
int *nframes, real **time,
int isize, atom_id index[],
real **trans_frac,
real **aver_angle,
real *dih[],
const output_env_t oenv)
{
t_pbc *pbc;
t_trxstatus *status;
int i, angind, natoms, total, teller;
int nangles, n_alloc;
real t, fraction, pifac, aa, angle;
real *angles[2];
matrix box;
rvec *x;
int cur = 0;
#define prev (1-cur)
snew(pbc, 1);
natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);
if (bAngles)
{
nangles = isize/3;
pifac = M_PI;
}
else
{
nangles = isize/4;
pifac = 2.0*M_PI;
}
snew(angles[cur], nangles);
snew(angles[prev], nangles);
/* Start the loop over frames */
total = 0;
teller = 0;
n_alloc = 0;
*time = NULL;
*trans_frac = NULL;
*aver_angle = NULL;
do
{
if (teller >= n_alloc)
{
n_alloc += 100;
if (bSaveAll)
{
for (i = 0; (i < nangles); i++)
{
srenew(dih[i], n_alloc);
}
}
srenew(*time, n_alloc);
srenew(*trans_frac, n_alloc);
srenew(*aver_angle, n_alloc);
}
(*time)[teller] = t;
if (pbc)
{
set_pbc(pbc, -1, box);
}
if (bAngles)
{
calc_angles(pbc, isize, index, angles[cur], x);
}
else
{
calc_dihs(pbc, isize, index, angles[cur], x);
/* Trans fraction */
fraction = calc_fraction(angles[cur], nangles);
(*trans_frac)[teller] = fraction;
/* Change Ryckaert-Bellemans dihedrals to polymer convention
* Modified 990913 by Erik:
* We actually shouldn't change the convention, since it's
* calculated from polymer above, but we change the intervall
* from [-180,180] to [0,360].
*/
if (bRb)
{
for (i = 0; (i < nangles); i++)
{
if (angles[cur][i] <= 0.0)
{
angles[cur][i] += 2*M_PI;
}
}
}
/* Periodicity in dihedral space... */
if (bPBC)
//.........这里部分代码省略.........
开发者ID:exianshine，项目名称:gromacs，代码行数:101，代码来源:anadih.c
示例18: gmx_dyecoupl
//.........这里部分代码省略.........
snew(khist, histbins);
}
do
{
clear_rvec(donvec);
clear_rvec(accvec);
clear_rvec(donpos);
clear_rvec(accpos);
for (i = 0; i < ndon / 2; i++)
{
rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
}
for (i = 0; i < nacc / 2; i++)
{
rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
}
unitv(donvec, donvec);
unitv(accvec, accvec);
svmul(1.0 / ndon, donpos, donpos);
svmul(1.0 / nacc, accpos, accpos);
if (bPBCdist)
{
set_pbc(pbc, ePBC, fr.box);
pbc_dx(pbc, donpos, accpos, dist);
}
else
{
rvec_sub(donpos, accpos, dist);
}
unitv(dist, distnorm);
R = norm(dist);
kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
kappa2 *= kappa2;
if (R0 > 0)
{
Rfrac = R/R0;
insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
insteffs += insteff;
if (bInstEffout)
{
fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
}
}
Rs += R;
kappa2s += kappa2;
rkcount++;
if (bRKout)
{
fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
}
开发者ID:rmcgibbo，项目名称:gromacs，代码行数:67，代码来源:gmx_dyecoupl.cpp
示例19: orient
void orient(int natom,rvec *x,rvec *v, rvec angle,matrix box)
{
real longest,rij,rzi;
int i,j,m,max_i=0,max_j=0;
rvec origin;
int temp;
real alfa=0,beta=0,gamma=0;
t_pbc pbc;
set_pbc(&pbc,-1,box);
/*first i am going to look for the longest atom-atom distance*/
longest=dist2(&pbc,x[0],x[1]);
i=0;
j=1;
for (i=0;(i<natom);i++) {
for (j=0;(j<natom);j++) {
rij=dist2(&pbc,x[i],x[j]);
if (rij>longest) {
max_i=i;
max_j=j;
longest=rij;
}
}
}
/* first check if x[max_i]<x[max_j] else swap*/
if (x[max_i][2]>x[max_j][2]) {
temp=max_i;
max_i=max_j;
max_j=temp;
}
/*set the origin to x[i]*/
for(m=0;(m<DIM);m++)
origin[m]=x[max_i][m];
for(i=0;(i<natom);i++)
for(m=0;(m<DIM);m++)
x[i][m]-=origin[m];
/* calculate the rotation angles alfa(x_axis) and beta(y_axis)
* the rotation angles must be calculated clockwise looking
* along the rotation axis to the origin*
* alfa (x-axis)
*/
alfa=atan(x[max_j][ZZ]/x[max_j][YY])-M_PI_2;
beta=M_PI_2-atan(x[max_j][ZZ]/x[max_j][XX]);
rotate_conf(natom,x,v,alfa,beta,gamma);
/* now search the longest distance for rotation along the z_axis */
longest=distance_to_z(x[0]);
max_i=0;
for (i=1;(i<natom);i++) {
rzi=distance_to_z(x[i]);
if (rzi>longest) {
longest = rzi;
max_i=i;
}
}
gamma=atan(x[max_i][YY]/x[max_i][XX])-M_PI_2;
rotate_conf(natom,x,v,0,0,gamma);
angle[0]=alfa;
angle[1]=beta;
angle[2]=gamma;
} /*orient()*/
开发者ID:TTarenzi，项目名称:MMCG-HAdResS，代码行数:64，代码来源:gbutil.c
示例20: do_force_lowlevel
void do_force_lowlevel(t_forcerec *fr, t_inputrec *ir,
t_idef *idef, t_commrec *cr,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
t_mdatoms *md,
rvec x[], history_t *hist,
rvec f[],
rvec f_longrange[],
gmx_enerdata_t *enerd,
t_fcdata *fcd,
gmx_localtop_t *top,
gmx_genborn_t *born,
gmx_bool bBornRadii,
matrix box,
t_lambda *fepvals,
real *lambda,
t_graph *graph,
t_blocka *excl,
rvec mu_tot[],
int flags,
float *cycles_pme)
{
int i, j;
int donb_flags;
gmx_bool bSB;
int pme_flags;
matrix boxs;
rvec box_size;
t_pbc pbc;
real dvdl_dum[efptNR], dvdl_nb[efptNR];
#ifdef GMX_MPI
double t0 = 0.0, t1, t2, t3; /* time measurement for coarse load balancing */
#endif
set_pbc(&pbc, fr->ePBC, box);
/* reset free energy components */
for (i = 0; i < efptNR; i++)
{
dvdl_nb[i] = 0;
dvdl_dum[i] = 0;
}
/* Reset box */
for (i = 0; (i < DIM); i++)
{
box_size[i] = box[i][i];
}
debug_gmx();
/* do QMMM first if requested */
if (fr->bQMMM)
{
enerd->term[F_EQM] = calculate_QMMM(cr, x, f, fr);
}
/* Call the short range functions all in one go. */
#ifdef GMX_MPI
/*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
#define TAKETIME FALSE
if (TAKETIME)
{
MPI_Barrier(cr->mpi_comm_mygroup);
t0 = MPI_Wtime();
}
#endif
if (ir->nwall)
{
/* foreign lambda component for walls */
real dvdl_walls = do_walls(ir, fr, box, md, x, f, lambda[efptVDW],
enerd->grpp.ener[egLJSR], nrnb);
enerd->dvdl_lin[efptVDW] += dvdl_walls;
}
/* If doing GB, reset dvda and calculate the Born radii */
if (ir->implicit_solvent)
{
wallcycle_sub_start(wcycle, ewcsNONBONDED);
for (i = 0; i < born->nr; i++)
{
fr->dvda[i] = 0;
}
if (bBornRadii)
{
calc_gb_rad(cr, fr, ir, top, x, &(fr->gblist), born, md, nrnb);
}
wallcycle_sub_stop(wcycle, ewcsNONBONDED);
}
where();
/* We only do non-bonded calculation with group scheme here, the verlet
* calls are done from do_force_cutsVERLET(). */
if (fr->cutoff_scheme == ecutsGROUP && (flags & GMX_FORCE_NONBONDED))
{
//.........这里部分代码省略.........
开发者ID:zlmturnout，项目名称:gromacs，代码行数:101，代码来源:force.cpp
注：本文中的set_pbc函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。 |
客服电话
APP下载
官方微信
请发表评论