float SculptIterateObject()

in layer2/Sculpt.cpp [1642:2426]


float SculptIterateObject(CSculpt * I, ObjectMolecule * obj,
                          int state, int n_cycle, float *center)
{
  PyMOLGlobals *G = I->G;
  CShaker *shk;
  int a0, a1, a2, a3, b0, b3;
  int aa;
  CoordSet *cs;
  float *disp = NULL;
  float *v, *v0, *v1, *v2, *v3;
  float diff[3], len;
  int *atm2idx = NULL;
  int *cnt = NULL;
  int *i;
  int hash;
  int nb_next;
  int h, k, l;
  int offset;
  float cutoff, vdw_cutoff;
  int ex;
  int eval_flag;
  int mask;
  float wt;
  float vdw;
  float vdw14;
  float vdw_wt;
  float vdw_wt14;
  float bond_wt;
  float angl_wt;
  float pyra_wt, pyra_inv_wt;
  float plan_wt;
  float line_wt;
  float tors_wt;
  float tors_tole;
  int active_flag = false;
  float hb_overlap, hb_overlap_base;
  int *active, n_active;
  int *exclude;
  AtomInfoType *ai0, *ai1;
  double task_time;
  float vdw_magnify, vdw_magnified = 1.0F;
  int nb_skip, nb_skip_count;
  float total_strain = 0.0F, strain;
  int total_count = 1;
  CGO *cgo = NULL;
  float good_color[3] = { 0.2, 1.0, 0.2 };
  float bad_color[3] = { 1.0, 0.2, 0.2 };
  int vdw_vis_mode;
  float vdw_vis_min = 0.0F, vdw_vis_mid = 0.0F, vdw_vis_max = 0.0F;
  float tri_sc, tri_wt;
  float min_sc, min_wt;
  float max_sc = 1.025F, max_wt = 0.75F;
  float *cs_coord;
  float solvent_radius;
  float avd_wt, avd_gp, avd_rg;
  int avd_ex;

  PRINTFD(G, FB_Sculpt)
    " SculptIterateObject-Debug: entered state=%d n_cycle=%d\n", state, n_cycle ENDFD;
  if(!n_cycle)
    n_cycle = -1;

  if((state < obj->NCSet) && obj->CSet[state] && n_cycle) {

    disp = Alloc(float, 3 * obj->NAtom);
    atm2idx = Alloc(int, obj->NAtom);
    cnt = Alloc(int, obj->NAtom);
    active = Alloc(int, obj->NAtom);
    exclude = Calloc(int, obj->NAtom);
    shk = I->Shaker;

    PRINTFD(G, FB_Sculpt)
      " SIO-Debug: NDistCon %d\n", shk->NDistCon ENDFD;

    cs = obj->CSet[state];
    cs_coord = cs->Coord;

    vdw = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_scale);
    vdw14 = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_scale14);
    vdw_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_weight);
    vdw_wt14 =
      SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_weight14);
    bond_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_bond_weight);
    angl_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_angl_weight);
    pyra_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_pyra_weight);
    pyra_inv_wt =
      SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_pyra_inv_weight);
    plan_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_plan_weight);
    line_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_line_weight);
    tri_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_tri_weight);
    tri_sc = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_tri_scale);

    min_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_min_weight);
    min_sc = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_min_scale);
    max_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_max_weight);
    max_sc = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_max_scale);

    mask = SettingGet_i(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_field_mask);
    hb_overlap =
      SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_hb_overlap);
    hb_overlap_base =
      SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_hb_overlap_base);
    tors_tole =
      SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_tors_tolerance);
    tors_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_tors_weight);
    vdw_vis_mode =
      SettingGet_i(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_vis_mode);
    solvent_radius =
      SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_solvent_radius);

    avd_wt = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_avd_weight);
    avd_gp = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_avd_gap);
    avd_rg = SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_avd_range);
    avd_ex = SettingGet_i(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_avd_excl);
    if(avd_gp < 0.0F)
      avd_gp = 1.5F * solvent_radius;
    if(avd_rg < 0.0F)
      avd_rg = solvent_radius;

    if(vdw_vis_mode) {
      vdw_vis_min =
        SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_vis_min);
      vdw_vis_mid =
        SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_vis_mid);
      vdw_vis_max =
        SettingGet_f(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_vdw_vis_max);

      if(!cs->SculptCGO)
        cs->SculptCGO = CGONew(G);
      else
        CGOReset(cs->SculptCGO);
    } else if(cs->SculptCGO) {
      CGOReset(cs->SculptCGO);
    }
    cgo = cs->SculptCGO;

    nb_skip = SettingGet_i(G, cs->Setting, obj->Obj.Setting, cSetting_sculpt_nb_interval);
    if(nb_skip > n_cycle)
      nb_skip = n_cycle;
    if(nb_skip < 0)
      nb_skip = 0;

    n_active = 0;
    ai0 = obj->AtomInfo;
    {
      int a;
      for(a = 0; a < obj->NAtom; a++) {
        if(ai0->flags & cAtomFlag_exclude) {
          exclude[a] = true;
          a1 = -1;
        } else {
          if(obj->DiscreteFlag) {
            if(cs == obj->DiscreteCSet[a]) {
              a1 = obj->DiscreteAtmToIdx[a];
            } else {
              a1 = -1;
            }
          } else {
            a1 = cs->AtmToIdx[a];
          }
        }
        if(a1 >= 0) {
          active_flag = true;
          active[n_active] = a;
          n_active++;
        }
        atm2idx[a] = a1;
        ai0++;
      }
    }

    if(active_flag) {

      /* first, create coordinate -> vertex mapping */
      /* and count number of constraints */

      task_time = UtilGetSeconds(G);
      vdw_magnify = 1.0F;
      nb_skip_count = 0;

      if(center) {
        int *a_ptr = active;
        int a;
        for(aa = 0; aa < n_active; aa++) {
          a = *(a_ptr++);
          {
            AtomInfoType *ai = obj->AtomInfo + a;
            if((ai->protekted != 1) && !(ai->flags & cAtomFlag_fix)) {
              v2 = cs_coord + 3 * atm2idx[a];
              center[4] += *(v2);
              center[5] += *(v2 + 1);
              center[6] += *(v2 + 2);
              center[7] += 1.0F;
            }
          }
        }
      }

      while(n_cycle--) {

        total_strain = 0.0F;
        total_count = 0;
        /* initialize displacements to zero */

        v = disp;
        i = cnt;
        for(aa = 0; aa < n_active; aa++) {
          int a = active[aa];
          v = disp + a * 3;
          cnt[a] = 0;
          *(v) = 0.0F;
          *(v + 1) = 0.0F;
          *(v + 2) = 0.0F;
        }

        /* apply distance constraints */

        {
          ShakerDistCon *sdc = shk->DistCon;
          int a,ndc = shk->NDistCon;
          for(a = 0; a < ndc; a++) {
            int sdc_type = sdc->type;
            int b1 = sdc->at0;
            int b2 = sdc->at1;

            switch (sdc_type) {
            case cShakerDistBond:
              eval_flag = cSculptBond & mask;
              wt = bond_wt;
              break;
            case cShakerDistAngle:
              eval_flag = cSculptAngl & mask;
              wt = angl_wt;
              break;
            case cShakerDistLimit:
              eval_flag = cSculptTri & mask;
              wt = tri_wt;
              break;
            case cShakerDistMinim:
              eval_flag = cSculptMin & mask;
              wt = min_wt * sdc->weight;
              break;
            case cShakerDistMaxim:
              eval_flag = cSculptMax & mask;
              wt = max_wt * sdc->weight;
              break;
            default:
              eval_flag = false;
              wt = 0.0F;
              break;
            }

            if(eval_flag && !(exclude[b1] || exclude[b2])) {
              a1 = atm2idx[b1]; /* coordinate set indices */
              a2 = atm2idx[b2];
              if((a1 >= 0) && (a2 >= 0)) {
                v1 = cs_coord + 3 * a1;
                v2 = cs_coord + 3 * a2;
                switch (sdc_type) {
                case cShakerDistLimit:
                  strain =
                    ShakerDoDistLimit(sdc->targ * tri_sc, v1, v2, disp + b1 * 3,
                                      disp + b2 * 3, wt);
                  if(strain > 0.0F) {
                    cnt[b1]++;
                    cnt[b2]++;
                    total_strain += strain;
                    total_count++;
                  }
                  break;
                case cShakerDistMaxim:
                  strain =
                    ShakerDoDistLimit(sdc->targ * max_sc, v1, v2, disp + b1 * 3,
                                      disp + b2 * 3, wt);
                  if(strain > 0.0F) {
                    cnt[b1]++;
                    cnt[b2]++;
                    total_strain += strain;
                    total_count++;
                  }
                  break;
                case cShakerDistMinim:
                  strain =
                    ShakerDoDistMinim(sdc->targ * min_sc, v1, v2, disp + b1 * 3,
                                      disp + b2 * 3, wt);
                  if(strain > 0.0F) {
                    cnt[b1]++;
                    cnt[b2]++;
                    total_strain += strain;
                    total_count++;
                  }
                  break;
                default:
                  total_strain +=
                    ShakerDoDist(sdc->targ, v1, v2, disp + b1 * 3, disp + b2 * 3, wt);
                  cnt[b1]++;
                  cnt[b2]++;
                  total_count++;
                }
              }
            }
            sdc++;
          }
        }
        /* apply line constraints */

        if(cSculptLine & mask) {
          ShakerLineCon *slc = shk->LineCon;
          int nlc = shk->NLineCon;
          int a,b1,b2;
          for(a = 0; a < nlc; a++) {
            b0 = slc->at0;
            b1 = slc->at1;
            b2 = slc->at2;
            a0 = atm2idx[b0];   /* coordinate set indices */
            a1 = atm2idx[b1];
            a2 = atm2idx[b2];

            if((a0 >= 0) && (a1 >= 0) && (a2 >= 0)
               && !(exclude[b0] || exclude[b1] || exclude[b2])) {
              cnt[b0]++;
              cnt[b1]++;
              cnt[b2]++;
              v0 = cs_coord + 3 * a0;
              v1 = cs_coord + 3 * a1;
              v2 = cs_coord + 3 * a2;
              total_strain +=
                ShakerDoLine(v0, v1, v2, disp + b0 * 3, disp + b1 * 3, disp + b2 * 3,
                             line_wt);
              total_count++;
            }
            slc++;
          }
        }

        /* apply pyramid constraints */

        if(cSculptPyra & mask) {
          ShakerPyraCon *spc = shk->PyraCon;
          int npc = shk->NPyraCon;
          int a,b1,b2;
          for(a = 0; a < npc; a++) {

            b0 = spc->at0;
            b1 = spc->at1;
            b2 = spc->at2;
            b3 = spc->at3;
            a0 = atm2idx[b0];
            a1 = atm2idx[b1];
            a2 = atm2idx[b2];
            a3 = atm2idx[b3];

            if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)
               && !(exclude[b0] || exclude[b1] || exclude[b2] || exclude[b3])) {
              v0 = cs_coord + 3 * a0;
              v1 = cs_coord + 3 * a1;
              v2 = cs_coord + 3 * a2;
              v3 = cs_coord + 3 * a3;
              total_strain += ShakerDoPyra(spc->targ1,
                                           spc->targ2,
                                           v0, v1, v2, v3,
                                           disp + b0 * 3,
                                           disp + b1 * 3,
                                           disp + b2 * 3,
                                           disp + b3 * 3, pyra_wt, pyra_inv_wt);
              total_count++;

              cnt[b0]++;
              cnt[b1]++;
              cnt[b2]++;
              cnt[b3]++;
            }
            spc++;
          }
        }

        if(cSculptPlan & mask) {
          ShakerPlanCon *snc = shk->PlanCon;
          int npc = shk->NPlanCon;
          int a,b1,b2;
          /* apply planarity constraints */

          snc = shk->PlanCon;
          for(a = 0; a < npc; a++) {

            b0 = snc->at0;
            b1 = snc->at1;
            b2 = snc->at2;
            b3 = snc->at3;
            a0 = atm2idx[b0];
            a1 = atm2idx[b1];
            a2 = atm2idx[b2];
            a3 = atm2idx[b3];

            if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)
               && !(exclude[b0] || exclude[b1] || exclude[b2] || exclude[b3])) {
              v0 = cs_coord + 3 * a0;
              v1 = cs_coord + 3 * a1;
              v2 = cs_coord + 3 * a2;
              v3 = cs_coord + 3 * a3;
              total_strain += ShakerDoPlan(v0, v1, v2, v3,
                                           disp + b0 * 3,
                                           disp + b1 * 3,
                                           disp + b2 * 3,
                                           disp + b3 * 3,
                                           snc->target, snc->fixed, plan_wt);
              total_count++;
              cnt[b0]++;
              cnt[b1]++;
              cnt[b2]++;
              cnt[b3]++;
            }

            snc++;
          }
        }

        /* apply torsion constraints */

        if(cSculptTors & mask) {
          ShakerTorsCon *stc = shk->TorsCon;
          int ntc = shk->NTorsCon;
          int a,b1,b2;
          /* apply planarity constraints */

          for(a = 0; a < ntc; a++) {

            b0 = stc->at0;
            b1 = stc->at1;
            b2 = stc->at2;
            b3 = stc->at3;
            a0 = atm2idx[b0];
            a1 = atm2idx[b1];
            a2 = atm2idx[b2];
            a3 = atm2idx[b3];

            if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)
               && !(exclude[b0] || exclude[b1] || exclude[b2] || exclude[b3])) {
              v0 = cs_coord + 3 * a0;
              v1 = cs_coord + 3 * a1;
              v2 = cs_coord + 3 * a2;
              v3 = cs_coord + 3 * a3;
              total_strain += ShakerDoTors(stc->type,
                                           v0, v1, v2, v3,
                                           disp + b0 * 3,
                                           disp + b1 * 3,
                                           disp + b2 * 3,
                                           disp + b3 * 3, tors_tole, tors_wt);
              total_count++;
              cnt[b0]++;
              cnt[b1]++;
              cnt[b2]++;
              cnt[b3]++;
            }
            stc++;
          }
        }
        /* apply nonbonded interactions */

        if((n_cycle > 0) && (nb_skip_count > 0)) {
          /*skip and then weight extra */
          nb_skip_count--;
          vdw_magnify += 1.0F;
        } else {
          int nb_off0, nb_off1;
          int v0i, v1i, v2i;
          int x0i;
          int don_b0;
          int acc_b0;
          int b1;
          vdw_magnified = vdw_magnify;
          vdw_magnify = 1.0F;

          nb_skip_count = nb_skip;
          if((cSculptVDW | cSculptVDW14 | cSculptAvoid) & mask) {
            /* compute non-bonded interations */

            /* construct nonbonded hash */

            nb_next = 1;
            for(aa = 0; aa < n_active; aa++) {
              b0 = active[aa];
              a0 = atm2idx[b0];
              VLACheck(I->NBList, int, nb_next + 2);
              v0 = cs_coord + 3 * a0;
              hash = nb_hash(v0);
              i = I->NBList + nb_next;
              *(i++) = *(I->NBHash + hash);
              *(i++) = hash;
              *(i++) = b0;
              *(I->NBHash + hash) = nb_next;
              nb_next += 3;
            }

            /* find neighbors for each atom */
            if((cSculptVDW | cSculptVDW14) & mask) {
              for(aa = 0; aa < n_active; aa++) {
                b0 = active[aa];
                a0 = atm2idx[b0];
                ai0 = obj->AtomInfo + b0;
                v0 = cs_coord + 3 * a0;
                don_b0 = I->Don[b0];
                acc_b0 = I->Acc[b0];
                v0i = (int) (*v0);
                v1i = (int) (*(v0 + 1));
                v2i = (int) (*(v0 + 2));
                x0i = ex_hash_i0(b0);
                for(h = -4; h < 5; h += 4) {
                  nb_off0 = nb_hash_off_i0(v0i, h);
                  for(k = -4; k < 5; k += 4) {
                    nb_off1 = nb_off0 | nb_hash_off_i1(v1i, k);
                    for(l = -4; l < 5; l += 4) {
                      /*  offset = *(I->NBHash+nb_hash_off(v0,h,k,l)); */
                      offset = *(I->NBHash + (nb_off1 | nb_hash_off_i2(v2i, l)));
                      while(offset) {
                        i = I->NBList + offset;
                        b1 = *(i + 2);
                        if(b1 > b0) {
                          /* determine exclusion (if any) */
                          {
                            int xoffset;
                            int *I_EXList = I->EXList;
                            int ex1;
                            int *j;
                            xoffset = *(I->EXHash + (x0i | ex_hash_i1(b1)));
                            ex = 10;
                            while(xoffset) {
                              xoffset = (*(j = I_EXList + xoffset));
                              if((*(j + 1) == b0) && (*(j + 2) == b1)) {
                                ex1 = *(j + 3);
                                if(ex1 < ex) {
                                  ex = ex1;
                                }
                              }
                            }
                          }
                          if(ex > 3) {
                            ai1 = obj->AtomInfo + b1;
                            cutoff = ai0->vdw + ai1->vdw;

                            if(ex == 10) {      /* standard interaction -- no exclusion */
                              if(don_b0 && I->Acc[b1]) {        /* h-bond */
                                if(ai0->protons == cAN_H) {
                                  cutoff -= hb_overlap;
                                } else {
                                  cutoff -= hb_overlap_base;
                                }
                              } else if(acc_b0 && I->Don[b1]) { /* h-bond */
                                if(ai1->protons == cAN_H) {
                                  cutoff -= hb_overlap;
                                } else {
                                  cutoff -= hb_overlap_base;
                                }
                              }
                              if(cSculptVDW & mask) {
                                vdw_cutoff = cutoff * vdw;
                                wt = vdw_wt * vdw_magnified;
                                a1 = atm2idx[b1];
                                v1 = cs_coord + 3 * a1;
                                if(vdw_vis_mode && cgo && (n_cycle < 1)
                                   && ((!((ai0->protekted && ai1->protekted)
                                          || (ai0->flags & ai1->flags & cAtomFlag_fix))
                                       ) || (ai0->flags & cAtomFlag_study)
                                       || (ai1->flags & cAtomFlag_study))) {
                                  SculptCGOBump(v0, v1, ai0->vdw, ai1->vdw, cutoff,
                                                vdw_vis_min, vdw_vis_mid, vdw_vis_max,
                                                good_color, bad_color, vdw_vis_mode, cgo);
                                }
                                if(SculptCheckBump(v0, v1, diff, &len, vdw_cutoff))
                                  if(SculptDoBump(vdw_cutoff, len, diff,
                                                  disp + b0 * 3, disp + b1 * 3, wt,
                                                  &total_strain)) {
                                    cnt[b0]++;
                                    cnt[b1]++;
                                    total_count++;
                                  }
                              }
                            } else if(ex == 4) {        /* 1-4 interation */
                              cutoff *= vdw14;
                              wt = vdw_wt14 * vdw_magnified;

                              if(cSculptVDW14 & mask) {
                                a1 = atm2idx[b1];
                                v1 = cs_coord + 3 * a1;
                                if(SculptCheckBump(v0, v1, diff, &len, cutoff)) {
                                  if(SculptDoBump(cutoff, len, diff,
                                                  disp + b0 * 3, disp + b1 * 3, wt,
                                                  &total_strain)) {
                                    cnt[b0]++;
                                    cnt[b1]++;
                                    total_count++;
                                  }
                                }
                              }
                            } else if(ex == 5) {
                              /* do nothing */
                            }
                          }
                        }
                        offset = (*i);
                      }
                    }
                  }
                }
              }
            }

            if(cSculptAvoid & mask) {
              float target;
              float range = solvent_radius * 0.75;
              /* tweak nb distances to avoid
                 sitting in the surface
                 rendition danger zone for too
                 long (vdw1+vdw2+0.75*solvent) */
              for(aa = 0; aa < n_active; aa++) {
                b0 = active[aa];
                a0 = atm2idx[b0];
                ai0 = obj->AtomInfo + b0;
                v0 = cs_coord + 3 * a0;
                don_b0 = I->Don[b0];
                acc_b0 = I->Acc[b0];
                v0i = (int) (*v0);
                v1i = (int) (*(v0 + 1));
                v2i = (int) (*(v0 + 2));
                x0i = ex_hash_i0(b0);
                for(h = -8; h < 9; h += 4) {
                  nb_off0 = nb_hash_off_i0(v0i, h);
                  for(k = -8; k < 9; k += 4) {
                    nb_off1 = nb_off0 | nb_hash_off_i1(v1i, k);
                    for(l = -8; l < 9; l += 4) {
                      /*  offset = *(I->NBHash+nb_hash_off(v0,h,k,l)); */
                      offset = *(I->NBHash + (nb_off1 | nb_hash_off_i2(v2i, l)));
                      while(offset) {
                        i = I->NBList + offset;
                        b1 = *(i + 2);
                        if(b1 > b0) {
                          /* determine exclusion (if any) */
                          {
                            int xoffset;
                            int *I_EXList = I->EXList;
                            int ex1;
                            int *j;
                            xoffset = *(I->EXHash + (x0i | ex_hash_i1(b1)));
                            ex = 10;
                            while(xoffset) {
                              xoffset = (*(j = I_EXList + xoffset));
                              if((*(j + 1) == b0) && (*(j + 2) == b1)) {
                                ex1 = *(j + 3);
                                if(ex1 < ex) {
                                  ex = ex1;
                                }
                              }
                            }
                          }
                          if(ex > avd_ex) {     /* either non-covalent or extended chain */
                            ai1 = obj->AtomInfo + b1;
                            target = ai0->vdw + ai1->vdw + avd_gp;
                            a1 = atm2idx[b1];
                            v1 = cs_coord + 3 * a1;

                            if(SculptCheckAvoid(v0, v1, diff, &len, target, avd_rg)) {
                              if(SculptDoAvoid(target, range, len, diff,
                                               disp + b0 * 3, disp + b1 * 3, avd_wt,
                                               &total_strain)) {
                                cnt[b0]++;
                                cnt[b1]++;
                                total_count++;
                              }
                            }
                          }
                        }
                        offset = (*i);
                      }
                    }
                  }
                }
              }
            }

            /* clean up nonbonded hash */

            i = I->NBList + 2;
            while(nb_next > 1) {
              *(I->NBHash + *i) = 0;
              i += 3;
              nb_next -= 3;
            }
          }
        }
        /* average the displacements */

        if(n_cycle >= 0) {
          int cnt_a,a;
          float _1 = 1.0F;
          float inv_cnt;
          int *a_ptr = active;
          float *lookup_inverse = I->inverse;
          for(aa = 0; aa < n_active; aa++) {
            if((cnt_a = cnt[(a = *(a_ptr++))])) {
              AtomInfoType *ai = obj->AtomInfo + a;
              RefPosType *cs_refpos = cs->RefPos;
              int flags;
              if(!(ai->protekted || ((flags = ai->flags) & cAtomFlag_fix))) {
                v1 = disp + 3 * a;
                v2 = cs_coord + 3 * atm2idx[a];

                if((flags & cAtomFlag_restrain) && cs_refpos) {
                  RefPosType *rp = cs_refpos + atm2idx[a];
                  if(rp->specified) {
                    float *v3 = rp->coord;
                    cnt_a++;
                    v1[0] += v3[0] - v2[0];
                    v1[1] += v3[1] - v2[1];
                    v1[2] += v3[2] - v2[2];
                  }
                }

                if(!(cnt_a & 0xFFFFFF00))       /* don't divide -- too slow! */
                  inv_cnt = lookup_inverse[cnt_a];
                else
                  inv_cnt = _1 / cnt_a;

                *(v2) += (*(v1)) * inv_cnt;
                *(v2 + 1) += (*(v1 + 1)) * inv_cnt;
                *(v2 + 2) += (*(v1 + 2)) * inv_cnt;
              }
            }
          }
          cs->invalidateRep(cRepAll, cRepInvCoord);
        } else if(cgo) {
          SceneDirty(G);
        }
        if(n_cycle <= 0) {
          int *a_ptr = active;
          if(center)
            for(aa = 0; aa < n_active; aa++) {
              int a = *(a_ptr++);
              {
                AtomInfoType *ai = obj->AtomInfo + a;
                if((ai->protekted != 1) && !(ai->flags & cAtomFlag_fix)) {
                  v2 = cs_coord + 3 * atm2idx[a];
                  center[0] += *(v2);
                  center[1] += *(v2 + 1);
                  center[2] += *(v2 + 2);
                  center[3] += 1.0F;
                }
              }
            }
          break;
        }
      }

      task_time = UtilGetSeconds(G) - task_time;
      PRINTFB(G, FB_Sculpt, FB_Blather)
        " Sculpt: %2.5f seconds %8.3f %d %8.3f\n", task_time, total_strain, total_count,
        100 * total_strain / total_count ENDFB(G);

      if(total_count)
        total_strain = (1000 * total_strain) / total_count;
    }
    FreeP(exclude);
    FreeP(active);
    FreeP(cnt);
    FreeP(disp);
    FreeP(atm2idx);
    if(cgo) {
      CGOStop(cgo);
      {
        int est = CGOCheckComplex(cgo);
        if(est) {
          cs->SculptCGO = CGOSimplify(cgo, est);
          CGOFree(cgo);
          CGOFree(cs->SculptShaderCGO);
        }
      }
    }

    EditorDihedralInvalid(G, obj);
  }

  PRINTFD(G, FB_Sculpt)
    " SculptIterateObject-Debug: leaving...\n" ENDFD;

  return total_strain;
}