本文整理汇总了C++中rvec_inc函数的典型用法代码示例。如果您正苦于以下问题：C++ rvec_inc函数的具体用法？C++ rvec_inc怎么用？C++ rvec_inc使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了rvec_inc函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: calc_com2
/* calculates com of all atoms in x[], *index has their index numbers
to get the masses from atom[] */
real calc_com2(rvec x[],int gnx,atom_id *index,t_mdatoms *md,rvec com,
matrix box)
{
int i,ii,m;
real m0,tm;
clear_rvec(com);
tm=0;
for(i=0; (i<gnx); i++) {
ii=index[i];
m0=md->massT[ii];
tm+=m0;
for(m=0; (m<DIM); m++)
com[m]+=m0*x[i][m];
}
svmul(1/tm,com,com);
for(m=DIM-1; m>=0; m--) {
/* next two lines used to be commented out */
if (com[m] < 0 ) rvec_inc(com,box[m]);
if (com[m] > box[m][m]) rvec_dec(com,box[m]);
}
return tm;
}示例2: orient_princ
void orient_princ(t_atoms *atoms, int isize, atom_id *index,
int natoms, rvec x[], rvec *v, rvec d)
{
int i, m;
rvec xcm, prcomp;
matrix trans;
calc_xcm(x, isize, index, atoms->atom, xcm, FALSE);
for (i = 0; i < natoms; i++)
{
rvec_dec(x[i], xcm);
}
principal_comp(isize, index, atoms->atom, x, trans, prcomp);
if (d)
{
copy_rvec(prcomp, d);
}
/* Check whether this trans matrix mirrors the molecule */
if (det(trans) < 0)
{
for (m = 0; (m < DIM); m++)
{
trans[ZZ][m] = -trans[ZZ][m];
}
}
rotate_atoms(natoms, NULL, x, trans);
if (v)
{
rotate_atoms(natoms, NULL, v, trans);
}
for (i = 0; i < natoms; i++)
{
rvec_inc(x[i], xcm);
}
}示例3: gmx_pme_receive_f
void gmx_pme_receive_f(t_commrec *cr,
rvec f[], matrix vir_q, real *energy_q,
matrix vir_lj, real *energy_lj,
real *dvdlambda_q, real *dvdlambda_lj,
float *pme_cycles)
{
int natoms, i;
#ifdef GMX_PME_DELAYED_WAIT
/* Wait for the x request to finish */
gmx_pme_send_coeffs_coords_wait(cr->dd);
#endif
natoms = cr->dd->nat_home;
if (natoms > cr->dd->pme_recv_f_alloc)
{
cr->dd->pme_recv_f_alloc = over_alloc_dd(natoms);
srenew(cr->dd->pme_recv_f_buf, cr->dd->pme_recv_f_alloc);
}
#ifdef GMX_MPI
MPI_Recv(cr->dd->pme_recv_f_buf[0],
natoms*sizeof(rvec), MPI_BYTE,
cr->dd->pme_nodeid, 0, cr->mpi_comm_mysim,
MPI_STATUS_IGNORE);
#endif
for (i = 0; i < natoms; i++)
{
rvec_inc(f[i], cr->dd->pme_recv_f_buf[i]);
}
receive_virial_energy(cr, vir_q, energy_q, vir_lj, energy_lj, dvdlambda_q, dvdlambda_lj, pme_cycles);
}示例4: reduce_thread_forces
static void reduce_thread_forces(int n, rvec *f,
tensor vir,
real *Vcorr,
int efpt_ind, real *dvdl,
int nthreads, f_thread_t *f_t)
{
int t, i;
/* This reduction can run over any number of threads */
#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntBonded)) private(t) schedule(static)
for (i = 0; i < n; i++)
{
for (t = 1; t < nthreads; t++)
{
rvec_inc(f[i], f_t[t].f[i]);
}
}
for (t = 1; t < nthreads; t++)
{
*Vcorr += f_t[t].Vcorr;
*dvdl += f_t[t].dvdl[efpt_ind];
m_add(vir, f_t[t].vir, vir);
}
}示例5: calc_geom
real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
rvec max, gmx_bool bDiam)
{
real diam2, d;
char *grpnames;
int ii, i, j;
clear_rvec(geom_center);
diam2 = 0;
if (isize == 0)
{
clear_rvec(min);
clear_rvec(max);
}
else
{
if (index)
{
ii = index[0];
}
else
{
ii = 0;
}
for (j = 0; j < DIM; j++)
{
min[j] = max[j] = x[ii][j];
}
for (i = 0; i < isize; i++)
{
if (index)
{
ii = index[i];
}
else
{
ii = i;
}
rvec_inc(geom_center, x[ii]);
for (j = 0; j < DIM; j++)
{
if (x[ii][j] < min[j])
{
min[j] = x[ii][j];
}
if (x[ii][j] > max[j])
{
max[j] = x[ii][j];
}
}
if (bDiam)
{
if (index)
{
for (j = i + 1; j < isize; j++)
{
d = distance2(x[ii], x[index[j]]);
diam2 = max(d, diam2);
}
}
else
{
for (j = i + 1; j < isize; j++)
{
d = distance2(x[i], x[j]);
diam2 = max(d, diam2);
}
}
}
}
svmul(1.0 / isize, geom_center, geom_center);
}
return sqrt(diam2);
}示例6: list_trn
static void list_trn(char *fn)
{
static real mass[5] = { 15.9994, 1.008, 1.008, 0.0, 0.0 };
int i,j=0,m,fpread,fpwrite,nframe;
rvec *x,*v,*f,fmol[2],xcm[2],torque[j],dx;
real mmm,len;
matrix box;
t_trnheader trn;
gmx_bool bOK;
printf("Going to open %s\n",fn);
fpread = open_trn(fn,"r");
fpwrite = open_tpx(NULL,"w");
gmx_fio_setdebug(fpwrite,TRUE);
mmm=mass[0]+2*mass[1];
for(i=0; (i<5); i++)
mass[i] /= mmm;
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)) {
if (trn.x_size && trn.f_size) {
printf("There are %d atoms\n",trn.natoms);
for(j=0; (j<2); j++) {
clear_rvec(xcm[j]);
clear_rvec(fmol[j]);
clear_rvec(torque[j]);
for(i=5*j; (i<5*j+5); i++) {
rvec_inc(fmol[j],f[i]);
for(m=0; (m<DIM); m++)
xcm[j][m] += mass[i%5]*x[i][m];
}
for(i=5*j; (i<5*j+5); i++) {
rvec_dec(x[i],xcm[j]);
cprod(x[i],f[i],dx);
rvec_inc(torque[j],dx);
rvec_inc(x[i],xcm[j]);
}
}
pr_rvecs(stdout,0,"FMOL ",fmol,2);
pr_rvecs(stdout,0,"TORQUE",torque,2);
printf("Distance matrix Water1-Water2\n%5s","");
for(j=0; (j<5); j++)
printf(" %10s",nm[j]);
printf("\n");
for(j=0; (j<5); j++) {
printf("%5s",nm[j]);
for(i=5; (i<10); i++) {
rvec_sub(x[i],x[j],dx);
len = sqrt(iprod(dx,dx));
printf(" %10.7f",len);
}
printf("\n");
}
}
}
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);
}示例7: gmx_traj
//.........这里部分代码省略.........
if (bOT && fr.bV)
{
fprintf(outt, " %g", time);
for (i = 0; i < ngroups; i++)
{
fprintf(outt, sffmt, temp(fr.v, mass, isize[i], index[i]));
}
fprintf(outt, "\n");
}
if (bEKT && fr.bV)
{
fprintf(outekt, " %g", time);
for (i = 0; i < ngroups; i++)
{
fprintf(outekt, sffmt, ektrans(fr.v, mass, isize[i], index[i]));
}
fprintf(outekt, "\n");
}
if (bEKR && fr.bX && fr.bV)
{
fprintf(outekr, " %g", time);
for (i = 0; i < ngroups; i++)
{
fprintf(outekr, sffmt, ekrot(fr.x, fr.v, mass, isize[i], index[i]));
}
fprintf(outekr, "\n");
}
if ((bCV || bCF) && fr.bX &&
(ctime < 0 || (fr.time >= ctime*0.999999 &&
fr.time <= ctime*1.000001)))
{
for (i = 0; i < fr.natoms; i++)
{
rvec_inc(sumx[i], fr.x[i]);
}
nr_xfr++;
}
if (bCV && fr.bV)
{
for (i = 0; i < fr.natoms; i++)
{
rvec_inc(sumv[i], fr.v[i]);
}
nr_vfr++;
}
if (bCF && fr.bF)
{
for (i = 0; i < fr.natoms; i++)
{
rvec_inc(sumf[i], fr.f[i]);
}
nr_ffr++;
}
}
while (read_next_frame(oenv, status, &fr));
if (gpbc != NULL)
{
gmx_rmpbc_done(gpbc);
}
/* clean up a bit */
close_trj(status);
if (bOX)示例8: gmx_covar
//.........这里部分代码省略.........
}
/* Prepare reference frame */
if (bPBC) {
gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,box);
gmx_rmpbc(gpbc,atoms->nr,box,xref);
}
if (bFit)
reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);
snew(x,natoms);
snew(xav,natoms);
ndim=natoms*DIM;
if (sqrt(GMX_LARGE_INT_MAX)<ndim) {
gmx_fatal(FARGS,"Number of degrees of freedoms to large for matrix.\n");
}
snew(mat,ndim*ndim);
fprintf(stderr,"Calculating the average structure ...\n");
nframes0 = 0;
nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
if (nat != atoms->nr)
fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat);
do {
nframes0++;
/* calculate x: a fitted struture of the selected atoms */
if (bPBC)
gmx_rmpbc(gpbc,nat,box,xread);
if (bFit) {
reset_x(nfit,ifit,nat,NULL,xread,w_rls);
do_fit(nat,w_rls,xref,xread);
}
for (i=0; i<natoms; i++)
rvec_inc(xav[i],xread[index[i]]);
} while (read_next_x(oenv,status,&t,nat,xread,box));
close_trj(status);
inv_nframes = 1.0/nframes0;
for(i=0; i<natoms; i++)
for(d=0; d<DIM; d++) {
xav[i][d] *= inv_nframes;
xread[index[i]][d] = xav[i][d];
}
write_sto_conf_indexed(opt2fn("-av",NFILE,fnm),"Average structure",
atoms,xread,NULL,epbcNONE,zerobox,natoms,index);
sfree(xread);
fprintf(stderr,"Constructing covariance matrix (%dx%d) ...\n",(int)ndim,(int)ndim);
nframes=0;
nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
tstart = t;
do {
nframes++;
tend = t;
/* calculate x: a (fitted) structure of the selected atoms */
if (bPBC)
gmx_rmpbc(gpbc,nat,box,xread);
if (bFit) {
reset_x(nfit,ifit,nat,NULL,xread,w_rls);
do_fit(nat,w_rls,xref,xread);
}
if (bRef)
for (i=0; i<natoms; i++)
rvec_sub(xread[index[i]],xref[index[i]],x[i]);
else
for (i=0; i<natoms; i++)示例9: dd_move_f_specat
void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
rvec *f, rvec *fshift)
{
gmx_specatsend_t *spas;
rvec *vbuf;
int n, n0, n1, d, dim, dir, i;
ivec vis;
int is;
gmx_bool bPBC, bScrew;
n = spac->at_end;
for (d = dd->ndim-1; d >= 0; d--)
{
dim = dd->dim[d];
if (dd->nc[dim] > 2)
{
/* Pulse the grid forward and backward */
spas = spac->spas[d];
n0 = spas[0].nrecv;
n1 = spas[1].nrecv;
n -= n1 + n0;
vbuf = spac->vbuf;
/* Send and receive the coordinates */
dd_sendrecv2_rvec(dd, d,
f+n+n1, n0, vbuf, spas[0].nsend,
f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
for (dir = 0; dir < 2; dir++)
{
bPBC = ((dir == 0 && dd->ci[dim] == 0) ||
(dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
bScrew = (bPBC && dd->bScrewPBC && dim == XX);
spas = &spac->spas[d][dir];
/* Sum the buffer into the required forces */
if (!bPBC || (!bScrew && fshift == NULL))
{
for (i = 0; i < spas->nsend; i++)
{
rvec_inc(f[spas->a[i]], *vbuf);
vbuf++;
}
}
else
{
clear_ivec(vis);
vis[dim] = (dir == 0 ? 1 : -1);
is = IVEC2IS(vis);
if (!bScrew)
{
/* Sum and add to shift forces */
for (i = 0; i < spas->nsend; i++)
{
rvec_inc(f[spas->a[i]], *vbuf);
rvec_inc(fshift[is], *vbuf);
vbuf++;
}
}
else
{
/* Rotate the forces */
for (i = 0; i < spas->nsend; i++)
{
f[spas->a[i]][XX] += (*vbuf)[XX];
f[spas->a[i]][YY] -= (*vbuf)[YY];
f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
if (fshift)
{
rvec_inc(fshift[is], *vbuf);
}
vbuf++;
}
}
}
}
}
else
{
/* Two cells, so we only need to communicate one way */
spas = &spac->spas[d][0];
n -= spas->nrecv;
/* Send and receive the coordinates */
dd_sendrecv_rvec(dd, d, dddirForward,
f+n, spas->nrecv, spac->vbuf, spas->nsend);
/* Sum the buffer into the required forces */
if (dd->bScrewPBC && dim == XX &&
(dd->ci[dim] == 0 ||
dd->ci[dim] == dd->nc[dim]-1))
{
for (i = 0; i < spas->nsend; i++)
{
/* Rotate the force */
f[spas->a[i]][XX] += spac->vbuf[i][XX];
f[spas->a[i]][YY] -= spac->vbuf[i][YY];
f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
}
}
else
{
for (i = 0; i < spas->nsend; i++)
{
//.........这里部分代码省略.........
示例10: do_listed_vdw_q
//.........这里部分代码省略.........
*/
if(ivdw==2)
{
gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
}
count = 0;
gmx_nb_free_energy_kernel(icoul,
ivdw,
i1,
&i0,
j_index,
&i1,
&shift_f,
fr->shift_vec[0],
fshift[0],
&gid,
x14[0],
f14[0],
chargeA,
chargeB,
eps,
krf,
crf,
fr->ewaldcoeff,
egcoul,
typeA,
typeB,
ntype,
nbfp,
egnb,
tabscale,
tab,
lambda,
dvdlambda,
fr->sc_alpha,
fr->sc_power,
fr->sc_sigma6_def,
fr->sc_sigma6_min,
TRUE,
&outeriter,
&inneriter);
}
else
{
/* Not perturbed - call kernel 330 */
nb_kernel330
( &i1,
&i0,
j_index,
&i1,
&shift_f,
fr->shift_vec[0],
fshift[0],
&gid,
x14[0],
f14[0],
chargeA,
&eps,
&krf,
&crf,
egcoul,
typeA,
&ntype,
nbfp,
egnb,
&tabscale,
tab,
NULL,
NULL,
NULL,
NULL,
&nthreads,
&count,
(void *)&mtx,
&outeriter,
&inneriter,
NULL);
}
/* Add the forces */
rvec_inc(f[ai],f14[0]);
rvec_dec(f[aj],f14[0]);
if (pf_global->bInitialized)
pf_atom_add_bonded(pf_global, ai, aj, PF_INTER_NB14, f14[0]);
if (g)
{
/* Correct the shift forces using the graph */
ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
shift_vir = IVEC2IS(dt);
rvec_inc(fshift[shift_vir],f14[0]);
rvec_dec(fshift[CENTRAL],f14[0]);
}
/* flops: eNR_KERNEL_OUTER + eNR_KERNEL330 + 12 */
}
}
return 0.0;
}示例11: unitcell
void unitcell(rvec x[],rvec box,gmx_bool bYaw,real odist,real hdist)
{
#define cx 0.81649658
#define cy 0.47140452
#define cy2 0.94280904
#define cz 0.33333333
rvec xx[24] = {
{ 0, 0, 0 }, /* O1 */
{ 0, 0, 1 }, /* H relative to Oxygen */
{ cx, cy, -cz },
{ cx, cy, -cz }, /* O2 */
{ 0, 0, -1 }, /* H relative to Oxygen */
{ cx,-cy, +cz },
{ cx, cy+cy2, 0 }, /* O3 */
{ -cx, cy, -cz }, /* H relative to Oxygen */
{ 0, -cy2, -cz },
{ 0, 2*cy+cy2, -cz }, /* O4 */
{-cx,-cy, +cz }, /* H relative to Oxygen */
{ 0 , cy2, +cz },
{ 0, 0, 1 }, /* O5 */
{-cx, cy, +cz }, /* H relative to Oxygen */
{ 0 , -cy2, +cz },
{ cx, cy, 1+cz }, /* O6 */
{ -cx, -cy, -cz }, /* H relative to Oxygen */
{ 0, cy2, -cz },
{ cx, cy+cy2, 1 }, /* O7 */
{ 0, 0, -1 }, /* H relative to Oxygen */
{ cx, cy, +cz },
{ 0, 2*cy+cy2,1+cz }, /* O8 */
{ 0, 0, 1 }, /* H relative to Oxygen */
{ cx, -cy, -cz }
};
int i,iin,iout,j,m;
rvec tmp,t2,dip;
clear_rvec(dip);
for(i=0; (i<8); i++) {
iin = 3*i;
if (bYaw)
iout = 5*i;
else
iout = iin;
svmul(odist,xx[iin],x[iout]);
svmul(-0.82,x[iout],t2);
rvec_inc(dip,t2);
for(j=1; (j<=2); j++) {
svmul(hdist,xx[iin+j],tmp);
rvec_add(x[iout],tmp,x[iout+j]);
svmul(0.41,x[iout+j],t2);
rvec_inc(dip,t2);
}
if (bYaw)
for(m=0; (m<DIM); m++)
x[iout+3][m] = x[iout+4][m] =
(1-2*DCONS)*x[iout][m]+DCONS*(x[iout+1][m]+x[iout+2][m]);
}
box[XX] = 2*cx;
box[YY] = 2*(cy2+cy);
box[ZZ] = 2*(1+cz);
for(i=0; (i<DIM); i++)
box[i] *= odist;
printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
}示例12: virial
void virial(FILE *fp,gmx_bool bFull,int nmol,rvec x[],matrix box,real rcut,
gmx_bool bYaw,real q[],gmx_bool bLJ)
{
int i,j,im,jm,natmol,ik,jk,m,ninter;
rvec dx,f,ftot,dvir,vir,pres,xcmi,xcmj,*force;
real dx6,dx2,dx1,fscal,c6,c12,vcoul,v12,v6,vctot,v12tot,v6tot;
t_pbc pbc;
set_pbc(&pbc,box);
fprintf(fp,"%3s - %3s: %6s %6s %6s %6s %8s %8s %8s\n",
"ai","aj","dx","dy","dz","|d|","virx","viry","virz");
clear_rvec(ftot);
clear_rvec(vir);
ninter = 0;
vctot = 0;
v12tot = 0;
v6tot = 0;
natmol = bYaw ? 5 : 3;
snew(force,nmol*natmol);
for(i=0; (i<nmol); i++) {
im = natmol*i;
/* Center of geometry */
clear_rvec(xcmi);
for(ik=0; (ik<natmol); ik++)
rvec_inc(xcmi,x[im+ik]);
for(m=0; (m<DIM); m++)
xcmi[m] /= natmol;
for(j=i+1; (j<nmol); j++) {
jm = natmol*j;
/* Center of geometry */
clear_rvec(xcmj);
for(jk=0; (jk<natmol); jk++)
rvec_inc(xcmj,x[jm+jk]);
for(m=0; (m<DIM); m++)
xcmj[m] /= natmol;
/* First check COM-COM distance */
pbc_dx(&pbc,xcmi,xcmj,dx);
if (norm(dx) < rcut) {
ninter++;
/* Neirest neighbour molecules! */
clear_rvec(dvir);
for(ik=0; (ik<natmol); ik++) {
for(jk=0; (jk<natmol); jk++) {
pbc_dx(&pbc,x[im+ik],x[jm+jk],dx);
dx2 = iprod(dx,dx);
dx1 = sqrt(dx2);
vcoul = q[ik]*q[jk]*ONE_4PI_EPS0/dx1;
vctot += vcoul;
if (bLJ) {
if (bYaw) {
c6 = yaw_lj[ik][2*jk];
c12 = yaw_lj[ik][2*jk+1];
}
else {
c6 = spc_lj[ik][2*jk];
c12 = spc_lj[ik][2*jk+1];
}
dx6 = dx2*dx2*dx2;
v6 = c6/dx6;
v12 = c12/(dx6*dx6);
v6tot -= v6;
v12tot+= v12;
fscal = (vcoul+12*v12-6*v6)/dx2;
}
else
fscal = vcoul/dx2;
for(m=0; (m<DIM); m++) {
f[m] = dx[m]*fscal;
dvir[m] -= 0.5*dx[m]*f[m];
}
rvec_inc(force[ik+im],f);
rvec_dec(force[jk+jm],f);
/*if (bFull)
fprintf(fp,"%3s%4d-%3s%4d: %6.3f %6.3f %6.3f %6.3f"
" %8.3f %8.3f %8.3f\n",
watname[ik],im+ik,watname[jk],jm+jk,
dx[XX],dx[YY],dx[ZZ],norm(dx),
dvir[XX],dvir[YY],dvir[ZZ]);*/
}
}
if (bFull)
fprintf(fp,"%3s%4d-%3s%4d: "
" %8.3f %8.3f %8.3f\n",
"SOL",i,"SOL",j,dvir[XX],dvir[YY],dvir[ZZ]);
rvec_inc(vir,dvir);
}
}
}
fprintf(fp,"There were %d interactions between the %d molecules (%.2f %%)\n",
ninter,nmol,(real)ninter/(0.5*nmol*(nmol-1)));
fprintf(fp,"Vcoul: %10.4e V12: %10.4e V6: %10.4e Vtot: %10.4e (kJ/mol)\n",
vctot/nmol,v12tot/nmol,v6tot/nmol,(vctot+v12tot+v6tot)/nmol);
pr_rvec(fp,0,"vir ",vir,DIM,TRUE);
for(m=0; (m<DIM); m++)
pres[m] = -2*PRESFAC/(det(box))*vir[m];
//.........这里部分代码省略.........
示例13: random_h_coords
void random_h_coords(int natmol,int nmol,rvec x[],rvec box,
gmx_bool bYaw,real odist,real hdist)
{
#define cx 0.81649658
#define cy 0.47140452
#define cy2 0.94280904
#define cz 0.33333333
rvec xx[24] = {
{ 0, 0, 0 }, /* O1 */
{ 0, 0, 1 }, /* H relative to Oxygen */
{ cx, cy, -cz },
{ cx, cy, -cz }, /* O2 */
{ 0, 0, -1 }, /* H relative to Oxygen */
{ cx,-cy, +cz },
{ cx, cy+cy2, 0 }, /* O3 */
{ -cx, cy, -cz }, /* H relative to Oxygen */
{ 0, -cy2, -cz },
{ 0, 2*cy+cy2, -cz }, /* O4 */
{-cx,-cy, +cz }, /* H relative to Oxygen */
{ 0 , cy2, +cz },
{ 0, 0, 1 }, /* O5 */
{-cx, cy, +cz }, /* H relative to Oxygen */
{ 0 , -cy2, +cz },
{ cx, cy, 1+cz }, /* O6 */
{ -cx, -cy, -cz }, /* H relative to Oxygen */
{ 0, cy2, -cz },
{ cx, cy+cy2, 1 }, /* O7 */
{ 0, 0, -1 }, /* H relative to Oxygen */
{ cx, cy, +cz },
{ 0, 2*cy+cy2,1+cz }, /* O8 */
{ 0, 0, 1 }, /* H relative to Oxygen */
{ cx, -cy, -cz }
};
int i,iin,iout,j,m;
rvec tmp,t2,dip;
clear_rvec(dip);
for(i=0; (i<nmol); i++) {
iin = natmol*i;
iout = iin;
svmul(odist,x[iin],x[iout]);
svmul(-0.82,x[iout],t2);
rvec_inc(dip,t2);
for(j=1; (j<=2); j++) {
svmul(hdist,xx[3*(i % 8)+j],tmp);
rvec_add(x[iout],tmp,x[iout+j]);
svmul(0.41,x[iout+j],t2);
rvec_inc(dip,t2);
}
}
box[XX] = 2*cx;
box[YY] = 2*(cy2+cy);
box[ZZ] = 2*(1+cz);
for(i=0; (i<DIM); i++)
box[i] *= odist;
printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
}示例14: ewald_LRcorrection
//.........这里部分代码省略.........
if (calc_excl_corr)
{
i1 = excl->index[i];
i2 = excl->index[i+1];
/* Loop over excluded neighbours */
for (j = i1; (j < i2); j++)
{
k = AA[j];
/*
* First we must test whether k <> i, and then, because the
* exclusions are all listed twice i->k and k->i we must select
* just one of the two.
* As a minor optimization we only compute forces when the charges
* are non-zero.
*/
if (k > i)
{
qqA = qiA*chargeA[k];
if (qqA != 0.0)
{
rvec_sub(x[i], x[k], dx);
if (bMolPBC)
{
/* Cheap pbc_dx, assume excluded pairs are at short distance. */
for (m = DIM-1; (m >= 0); m--)
{
if (dx[m] > 0.5*box[m][m])
{
rvec_dec(dx, box[m]);
}
else if (dx[m] < -0.5*box[m][m])
{
rvec_inc(dx, box[m]);
}
}
}
dr2 = norm2(dx);
/* Distance between two excluded particles may be zero in the
* case of shells
*/
if (dr2 != 0)
{
rinv = gmx_invsqrt(dr2);
rinv2 = rinv*rinv;
dr = 1.0/rinv;
#ifdef TABLES
r1t = tabscale*dr;
n0 = r1t;
assert(n0 >= 3);
n1 = 12*n0;
eps = r1t-n0;
eps2 = eps*eps;
nnn = n1;
Y = VFtab[nnn];
F = VFtab[nnn+1];
Geps = eps*VFtab[nnn+2];
Heps2 = eps2*VFtab[nnn+3];
Fp = F+Geps+Heps2;
VV = Y+eps*Fp;
FF = Fp+Geps+2.0*Heps2;
vc = qqA*(rinv-VV);
fijC = qqA*FF;
Vexcl += vc;
fscal = vc*rinv2+fijC*tabscale*rinv;示例15: dielectric
//.........这里部分代码省略.........
}
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);
if (bINT)
{
srenew(djc, valloc);
}
srenew(v0, valloc);
if (bACF)
{
srenew(cacf, valloc);
}
}
if (time[nfr] <= bvit)
{
ii = nvfr;