/* A* ------------------------------------------------------------------- B* This file contains source code for the PyMOL computer program C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific. D* ------------------------------------------------------------------- E* It is unlawful to modify or remove this copyright notice. F* ------------------------------------------------------------------- G* Please see the accompanying LICENSE file for further information. H* ------------------------------------------------------------------- I* Additional authors of this source file include: -* -* -* Z* ------------------------------------------------------------------- */ #include"os_predef.h" #include"os_std.h" #include"os_gl.h" #include"Base.h" #include"Debug.h" #include"Parse.h" #include"OOMac.h" #include"Vector.h" #include"MemoryDebug.h" #include"Err.h" #include"Map.h" #include"Selector.h" #include"ObjectMolecule.h" #include"Ortho.h" #include"Util.h" #include"Vector.h" #include"Selector.h" #include"Matrix.h" #include"Scene.h" #include"P.h" #include"PConv.h" #include"Executive.h" #include"Setting.h" #include"Sphere.h" #include"main.h" #include"CGO.h" #include"Raw.h" #include"Editor.h" #include"Selector.h" #include"Sculpt.h" #define cMaxNegResi 100 #define ntrim ParseNTrim #define nextline ParseNextLine #define ncopy ParseNCopy #define nskip ParseNSkip void ObjectMoleculeRender(ObjectMolecule *I,int frame,CRay *ray,Pickable **pick,int pass); void ObjectMoleculeCylinders(ObjectMolecule *I); CoordSet *ObjectMoleculeMMDStr2CoordSet(char *buffer,AtomInfoType **atInfoPtr); int ObjectMoleculeDoesAtomNeighborSele(ObjectMolecule *I, int index, int sele); CoordSet *ObjectMoleculePMO2CoordSet(CRaw *pmo,AtomInfoType **atInfoPtr,int *restart); void ObjectMoleculeAppendAtoms(ObjectMolecule *I,AtomInfoType *atInfo,CoordSet *cset); void ObjectMoleculeUpdate(ObjectMolecule *I); int ObjectMoleculeGetNFrames(ObjectMolecule *I); void ObjectMoleculeDescribeElement(ObjectMolecule *I,int index,char *buffer); void ObjectMoleculeSeleOp(ObjectMolecule *I,int sele,ObjectMoleculeOpRec *op); void ObjectMoleculeTransformTTTf(ObjectMolecule *I,float *ttt,int state); CoordSet *ObjectMoleculeChemPyModel2CoordSet(PyObject *model,AtomInfoType **atInfoPtr); int ObjectMoleculeGetAtomGeometry(ObjectMolecule *I,int state,int at); void ObjectMoleculeBracketResidue(ObjectMolecule *I,AtomInfoType *ai,int *st,int *nd); void ObjectMoleculeAddSeleHydrogens(ObjectMolecule *I,int sele); CoordSet *ObjectMoleculeXYZStr2CoordSet(char *buffer,AtomInfoType **atInfoPtr); CSetting **ObjectMoleculeGetSettingHandle(ObjectMolecule *I,int state); void ObjectMoleculeInferAmineGeomFromBonds(ObjectMolecule *I,int state); CoordSet *ObjectMoleculeTOPStr2CoordSet(char *buffer, AtomInfoType **atInfoPtr); ObjectMolecule *ObjectMoleculeReadTOPStr(ObjectMolecule *I,char *TOPStr,int frame,int discrete); void ObjectMoleculeInferHBondFromChem(ObjectMolecule *I); static char *ObjectMoleculeGetCaption(ObjectMolecule *I) { int state = ObjectGetCurrentState((CObject*)I,false); if((state>=0)&&(state<I->NCSet)) { CoordSet *cs = I->CSet[state]; if(cs) { return cs->Name; } } return NULL; } #define MAX_BOND_DIST 50 #if 0 static void dump_jxn(char *lab,char *q) { char cc[MAXLINELEN]; printf("\n%s %p\n",lab,q); q=q-150; q=nextline(q); ncopy(cc,q,100); printf("0[%s]\n",cc); q=nextline(q); ncopy(cc,q,100); printf("1[%s]\n",cc); q=nextline(q); ncopy(cc,q,100); printf("2[%s]\n",cc); q=nextline(q); ncopy(cc,q,100); printf("3[%s]\n",cc); } #endif static char *skip_fortran(int num,int per_line,char *p) { int a,b; b=0; for(a=0;a<num;a++) { if((++b)==per_line) { b=0; p=nextline(p); } } if(b) p=nextline(p); return(p); } void M4XAnnoInit(M4XAnnoType *m4x) { UtilZeroMem((char*)m4x,sizeof(M4XAnnoType)); } void M4XAnnoPurge(M4XAnnoType *m4x) { int c; if(m4x) { for(c=0;c<m4x->n_context;c++) { VLAFreeP(m4x->context[c].hbond); VLAFreeP(m4x->context[c].nbond); VLAFreeP(m4x->context[c].site); VLAFreeP(m4x->context[c].ligand); VLAFreeP(m4x->context[c].water); } if(m4x->align) { M4XAlignPurge(m4x->align); } VLAFreeP(m4x->context); } } void M4XAlignInit(M4XAlignType *align) { UtilZeroMem((char*)align,sizeof(M4XAlignType)); align->id_at_point = VLACalloc(int,100); align->fitness = VLAlloc(float,100); } void M4XAlignPurge(M4XAlignType *align) { VLAFreeP(align->id_at_point); VLAFreeP(align->fitness); FreeP(align); } void ObjectMoleculeOpRecInit(ObjectMoleculeOpRec *op) { UtilZeroMem((char*)op,sizeof(ObjectMoleculeOpRec)); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadTRJFile(ObjectMolecule *I,char *fname,int frame, int interval,int average,int start, int stop,int max,char *sele,int image, float *shift) { int ok=true; FILE *f; char *buffer,*p; char cc[MAXLINELEN]; int n_read; int to_go; int skip_first_line = true; int periodic=false; int angles=true; float f0,f1,f2,*fp; float box[3],angle[3]; float r_cent[3],r_trans[3]; int r_act,r_val,r_cnt; float *r_fp_start=NULL,*r_fp_stop=NULL; int a,b,c,i; int *to,at_i; int zoom_flag=false; int cnt=0; int n_avg=0; int icnt; int ncnt=0; int sele0 = SelectorIndexByName(sele); int *xref = NULL; float zerovector[3]={0.0,0.0,0.0}; CoordSet *cs = NULL; if(!shift) shift = zerovector; if(interval<1) interval=1; icnt=interval; #define BUFSIZE 4194304 #define GETTING_LOW 10000 f=fopen(fname,"rb"); if(!f) { ok=ErrMessage("ObjectMoleculeLoadTOPFile","Unable to open file!"); } else { if(!I->CSTmpl) { PRINTFB(FB_Errors,FB_ObjectMolecule) " ObjMolLoadTRJFile: Missing topology" ENDFB; return(I); } cs=CoordSetCopy(I->CSTmpl); if(sele0>=0) { /* build array of cross-references */ xref = Alloc(int,I->NAtom); c=0; for(a=0;a<I->NAtom;a++) { if(SelectorIsMember(I->AtomInfo[a].selEntry,sele0)) { xref[a]=c++; } else { xref[a]=-1; } } for(a=0;a<I->NAtom;a++) { /* now terminate the excluded references */ if(xref[a]<0) { cs->AtmToIdx[a]=-1; } } to=cs->IdxToAtm; c=0; for(a=0;a<cs->NIndex;a++) { /* now fix IdxToAtm, AtmToIdx, and remap xref to coordinate space */ at_i = cs->IdxToAtm[a]; if(cs->AtmToIdx[at_i]>=0) { *(to++)=at_i; cs->AtmToIdx[at_i]=c; xref[a]=c; c++; } else { xref[a]=-1; } } cs->NIndex=c; cs->IdxToAtm = Realloc(cs->IdxToAtm,int,cs->NIndex+1); VLASize(cs->Coord,float,cs->NIndex*3); } PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjMolLoadTRJFile: Loading from \"%s\".\n",fname ENDFB; buffer = (char*)mmalloc(BUFSIZE+1); /* 1 MB read buffer */ p = buffer; buffer[0]=0; n_read = 0; to_go=0; a = 0; b = 0; c = 0; f1=0.0; f2=0.0; while(1) { to_go = n_read-(p-buffer); if(to_go<GETTING_LOW) if(!feof(f)) { if(to_go) memcpy(buffer,p,to_go); n_read = fread(buffer+to_go,1,BUFSIZE-to_go,f); n_read = to_go + n_read; buffer[n_read]=0; p = buffer; if(skip_first_line) { p=nextline(p); skip_first_line=false; } to_go = n_read-(p-buffer); } if(!to_go) break; p=ncopy(cc,p,8); if((++b)==10) { b=0; p=nextline(p); } f0 = f1; f1 = f2; if(sscanf(cc,"%f",&f2)==1) { if((++c)==3) { c=0; if((cnt+1)>=start) { if(icnt<=1) { if(xref) { if(xref[a]>=0) fp=cs->Coord+3*xref[a]; else fp=NULL; } else { fp=cs->Coord+3*a; } if(fp) { if(n_avg) { *(fp++)+=f0; *(fp++)+=f1; *(fp++)+=f2; } else { *(fp++)=f0; *(fp++)=f1; *(fp++)=f2; } } } } if((++a)==I->NAtom) { cnt++; a=0; if(b) p=nextline(p); b=0; if(cs->PeriodicBoxType!=cCSet_NoPeriodicity) { /* read periodic box */ c=0; periodic=true; angles=true; p = ncopy(cc,p,8); if(sscanf(cc,"%f",&box[0])!=1) periodic=false; p = ncopy(cc,p,8); if(sscanf(cc,"%f",&box[1])!=1) periodic=false; p = ncopy(cc,p,8); if(sscanf(cc,"%f",&box[2])!=1) periodic=false; p = ncopy(cc,p,8); if(sscanf(cc,"%f",&angle[0])!=1) angles=false; p = ncopy(cc,p,8); if(sscanf(cc,"%f",&angle[1])!=1) angles=false; p = ncopy(cc,p,8); if(sscanf(cc,"%f",&angle[2])!=1) angles=false; if(periodic) { if(!cs->PeriodicBox) cs->PeriodicBox=CrystalNew(); cs->PeriodicBox->Dim[0] = box[0]; cs->PeriodicBox->Dim[1] = box[1]; cs->PeriodicBox->Dim[2] = box[2]; if(angles) { cs->PeriodicBox->Angle[0] = angle[0]; cs->PeriodicBox->Angle[1] = angle[1]; cs->PeriodicBox->Angle[2] = angle[2]; } CrystalUpdate(cs->PeriodicBox); /* CrystalDump(cs->PeriodicBox);*/ p=nextline(p); b=0; } if(cs->PeriodicBoxType==cCSet_Octahedral) periodic=false; /* can't handle this yet... */ } if((stop>0)&&(cnt>=stop)) break; if(cnt>=start) { icnt--; if(icnt>0) { PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: skipping set %d...\n",cnt ENDFB; } else { icnt=interval; n_avg++; } if(icnt==interval) { if(n_avg<average) { PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: averaging set %d...\n",cnt ENDFB; } else { /* compute average */ if(n_avg>1) { fp=cs->Coord; for(i=0;i<cs->NIndex;i++) { *(fp++)/=n_avg; *(fp++)/=n_avg; *(fp++)/=n_avg; } } if(periodic&&image) { /* Perform residue-based period image transformation */ i = 0; r_cnt = 0; r_act = 0; /* 0 unspec, 1=load, 2=image, 3=leave*/ r_val = -1; while(r_act!=3) { if(i>=cs->NIndex) { if(r_cnt) r_act = 2; else r_act = 3; } if(r_act==0) { /* start new residue */ r_cnt = 0; r_act = 1; /* now load */ } if(r_act==1) { if(i<cs->NIndex) { /* is there a coordinate for atom? */ if(xref) { if(xref[i]>=0) fp=cs->Coord+3*xref[i]; else fp=NULL; } else { fp=cs->Coord+3*i; } if(fp) { /* yes there is... */ if(r_cnt) { if(r_val!=I->AtomInfo[cs->IdxToAtm[i]].resv) { r_act=2; /* end of residue-> time to image */ } else { r_cnt++; r_cent[0]+=*(fp++); r_cent[1]+=*(fp++); r_cent[2]+=*(fp++); r_fp_stop = fp; /* stop here */ i++; } } else { r_val = I->AtomInfo[cs->IdxToAtm[i]].resv; r_cnt++; r_fp_start = fp; /* start here */ r_cent[0]=*(fp++); r_cent[1]=*(fp++); r_cent[2]=*(fp++); r_fp_stop = fp; /* stop here */ i++; } } else { i++; } } else { r_act=2; /* image */ } } if(r_act==2) { /* time to image */ if(r_cnt) { r_cent[0]/=r_cnt; r_cent[1]/=r_cnt; r_cent[2]/=r_cnt; transform33f3f(cs->PeriodicBox->RealToFrac,r_cent,r_cent); r_trans[0]=(float)fmod(1000.0+shift[0]+r_cent[0],1.0F); r_trans[1]=(float)fmod(1000.0+shift[1]+r_cent[1],1.0F); r_trans[2]=(float)fmod(1000.0+shift[2]+r_cent[2],1.0F); r_trans[0]-=r_cent[0]; r_trans[1]-=r_cent[1]; r_trans[2]-=r_cent[2]; transform33f3f(cs->PeriodicBox->FracToReal,r_trans,r_trans); fp=r_fp_start; while(fp<r_fp_stop) { *(fp++)+=r_trans[0]; *(fp++)+=r_trans[1]; *(fp++)+=r_trans[2]; } } r_act=0; /* reset */ r_cnt=0; } } } /* add new coord set */ if(cs->fInvalidateRep) cs->fInvalidateRep(cs,cRepAll,cRepInvRep); if(frame<0) frame=I->NCSet; if(!I->NCSet) { zoom_flag=true; } VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cs; ncnt++; if(average<2) { PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: read set %d into state %d...\n",cnt,frame+1 ENDFB; } else { PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: averaging set %d...\n",cnt ENDFB; PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: average loaded into state %d...\n",frame+1 ENDFB; } frame++; cs = CoordSetCopy(cs); n_avg=0; if((stop>0)&&(cnt>=stop)) break; if((max>0)&&(ncnt>=max)) break; } } } else { PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: skipping set %d...\n",cnt ENDFB; } } } } else { PRINTFB(FB_Errors,FB_ObjectMolecule) " ObjMolLoadTRJFile-Error: Failed to read an expected coordinate value.\n This trajectory does not match the loaded parameter/topology file.\n Likely cause: either the atom count or the periodic box settings\n are inconsistent between the two files.\n" ENDFB; break; } } FreeP(xref); mfree(buffer); } if(cs) cs->fFree(cs); SceneChanged(); SceneCountFrames(); if(zoom_flag) if(SettingGet(cSetting_auto_zoom)) { ExecutiveWindowZoom(I->Obj.Name,0.0,-1,0); /* auto zoom (all states) */ } return(I); } ObjectMolecule *ObjectMoleculeLoadRSTFile(ObjectMolecule *I,char *fname,int frame) { int ok=true; FILE *f; char *buffer,*p; char cc[MAXLINELEN]; float f0,f1,f2,*fp; int a,b,c; int zoom_flag=false; CoordSet *cs = NULL; int size; #define BUFSIZE 4194304 #define GETTING_LOW 10000 f=fopen(fname,"rb"); if(!f) ok=ErrMessage("ObjectMoleculeLoadRSTFile","Unable to open file!"); else { if(!I->CSTmpl) { PRINTFB(FB_Errors,FB_ObjectMolecule) " ObjMolLoadTRJFile: Missing topology" ENDFB; return(I); } cs=CoordSetCopy(I->CSTmpl); PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjMolLoadTRJFile: Loading from \"%s\".\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); p=nextline(p); p=nextline(p); a = 0; b = 0; c = 0; f1=0.0; f2=0.0; while(*p) { p=ncopy(cc,p,12); if((++b)==6) { b=0; p=nextline(p); } f0 = f1; f1 = f2; if(sscanf(cc,"%f",&f2)==1) { if((++c)==3) { c=0; fp=cs->Coord+3*a; *(fp++)=f0; *(fp++)=f1; *(fp++)=f2; if((++a)==I->NAtom) { a=0; if(b) p=nextline(p); b=0; /* add new coord set */ if(cs->fInvalidateRep) cs->fInvalidateRep(cs,cRepAll,cRepInvRep); if(frame<0) frame=I->NCSet; if(!I->NCSet) { zoom_flag=true; } VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cs; PRINTFB(FB_Details,FB_ObjectMolecule) " ObjectMolecule: read coordinates into state %d...\n",frame+1 ENDFB; cs = CoordSetCopy(cs); break; } } } else { PRINTFB(FB_Errors,FB_ObjectMolecule) " ObjMolLoadTRJFile: atom/coordinate mismatch.\n" ENDFB; break; } } mfree(buffer); } if(cs) cs->fFree(cs); SceneChanged(); SceneCountFrames(); if(zoom_flag) if(SettingGet(cSetting_auto_zoom)) { ExecutiveWindowZoom(I->Obj.Name,0.0,-1,0); /* auto zoom (all states) */ } return(I); } static char *findflag(char *p,char *flag,char *format) { char cc[MAXLINELEN]; char pat[MAXLINELEN] = "%FLAG "; int l; PRINTFD(FB_ObjectMolecule) " findflag: flag %s format %s\n",flag,format ENDFD; strcat(pat,flag); l=strlen(pat); while(*p) { p=ncopy(cc,p,l); if(WordMatch(cc,pat,true)<0) { p=nextline(p); break; } p=nextline(p); if(!*p) { PRINTFB(FB_ObjectMolecule,FB_Errors) " ObjectMolecule-Error: Unrecognized file format (can't find \"%s\").\n",pat ENDFB; } } strcpy(pat,"%FORMAT("); strcat(pat,format); strcat(pat,")"); l=strlen(pat); while(*p) { p=ncopy(cc,p,l); if(WordMatch(cc,pat,true)<0) { p=nextline(p); break; } p=nextline(p); if(!*p) { PRINTFB(FB_ObjectMolecule,FB_Errors) " ObjectMolecule-Error: Unrecognized file format (can't find \"%s\").\n",pat ENDFB; } } return(p); } #define nextline_top nextline /*========================================================================*/ CoordSet *ObjectMoleculeTOPStr2CoordSet(char *buffer, AtomInfoType **atInfoPtr) { char *p; int nAtom; int a,b,c,bi,last_i,at_i,aa,rc; float *coord = NULL; float *f; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL,*ai; BondType *bond=NULL,*bd; int nBond=0; int auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); int auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); int amber7 = false; WordType title; ResName *resn; char cc[MAXLINELEN]; int ok=true; int i0,i1,i2; /* trajectory parameters */ int NTYPES,NBONH,MBONA,NTHETH,MTHETA; int NPHIH,MPHIA,NHPARM,NPARM,NNB,NRES; int NBONA,NTHETA,NPHIA,NUMBND,NUMANG,NPTRA; int NATYP,NPHB,IFPERT,NBPER,NGPER,NDPER; int MBPER,MGPER,MDPER,IFBOX,NMXRS,IFCAP; int NEXTRA,IPOL=0; int wid,col; float BETA; float BOX1,BOX2,BOX3; cset = CoordSetNew(); p=buffer; nAtom=0; if(atInfoPtr) atInfo = *atInfoPtr; if(!atInfo) ErrFatal("TOPStr2CoordSet","need atom information record!"); /* failsafe for old version..*/ ncopy(cc,p,8); if(strcmp(cc,"%VERSION")==0) { amber7=true; PRINTFB(FB_ObjectMolecule,FB_Details) " ObjectMolecule: Attempting to read Amber7 topology file.\n" ENDFB; } else { PRINTFB(FB_ObjectMolecule,FB_Details) " ObjectMolecule: Assuming this is an Amber6 topology file.\n" ENDFB; } /* read title */ if(amber7) { p = findflag(buffer,"TITLE","20a4"); } p=ncopy(cc,p,20); title[0]=0; sscanf(cc,"%s",title); p=nextline_top(p); if(amber7) { p = findflag(buffer,"POINTERS","10I8"); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&nAtom)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NTYPES)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NBONH)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&MBONA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NTHETH)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&MTHETA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NPHIH)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&MPHIA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NHPARM)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NPARM)==1); p=nextline_top(p); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NNB)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NRES)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NBONA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NTHETA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NPHIA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NUMBND)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NUMANG)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NPTRA)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NATYP)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NPHB)==1); p=nextline_top(p); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&IFPERT)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NBPER)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NGPER)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NDPER)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&MBPER)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&MGPER)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&MDPER)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&IFBOX)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NMXRS)==1); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&IFCAP)==1); p=nextline_top(p); p=ncopy(cc,p,8); ok = ok && (sscanf(cc,"%d",&NEXTRA)==1); } else { p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&nAtom)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NTYPES)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NBONH)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&MBONA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NTHETH)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&MTHETA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NPHIH)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&MPHIA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NHPARM)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NPARM)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NNB)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NRES)==1); p=nextline_top(p); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NBONA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NTHETA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NPHIA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NUMBND)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NUMANG)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NPTRA)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NATYP)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NPHB)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&IFPERT)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NBPER)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NGPER)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NDPER)==1); p=nextline_top(p); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&MBPER)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&MGPER)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&MDPER)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&IFBOX)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&NMXRS)==1); p=ncopy(cc,p,6); ok = ok && (sscanf(cc,"%d",&IFCAP)==1); p=ncopy(cc,p,6); if(sscanf(cc,"%d",&NEXTRA)!=1) NEXTRA=0; } switch(IFBOX) { case 2: cset->PeriodicBoxType = cCSet_Octahedral; PRINTFB(FB_ObjectMolecule,FB_Details) " TOPStrToCoordSet: Warning: can't currently image a truncated octahedron...\n" ENDFB; break; case 1: cset->PeriodicBoxType = cCSet_Orthogonal; break; case 0: default: cset->PeriodicBoxType = cCSet_NoPeriodicity; break; } p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error reading counts lines"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read counts line nAtom %d NBONA %d NBONH %d\n", nAtom,NBONA,NBONH ENDFB; } if(ok) { VLACheck(atInfo,AtomInfoType,nAtom); if(amber7) { p = findflag(buffer,"ATOM_NAME","20a4"); } /* read atoms */ b=0; for(a=0;a<nAtom;a++) { p=ncopy(cc,p,4); ai=atInfo+a; if(!sscanf(cc,"%s",ai->name)) ai->name[0]=0; if((++b)==20) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error reading atom names"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read atom names.\n" ENDFB; } /* read charges */ if(amber7) { p = findflag(buffer,"CHARGE","5E16.8"); } b=0; for(a=0;a<nAtom;a++) { p=ncopy(cc,p,16); ai=atInfo+a; if(!sscanf(cc,"%f",&ai->partialCharge)) ok=false; else { ai->partialCharge/=18.2223F; /* convert to electron charge */ } if((++b)==5) { b=0; p=nextline_top(p); } } if(!ok) { ErrMessage("TOPStrToCoordSet","Error reading charges"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read charges.\n" ENDFB; } if(b) p=nextline_top(p); if(!amber7) { /* skip masses */ p=skip_fortran(nAtom,5,p); } /* read LJ atom types */ if(amber7) { p = findflag(buffer,"ATOM_TYPE_INDEX","10I8"); col=10; wid=8; } else { col=12; wid=6; } b=0; for(a=0;a<nAtom;a++) { p=ncopy(cc,p,wid); ai=atInfo+a; if(!sscanf(cc,"%d",&ai->customType)) ok=false; if((++b)==col) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error LJ atom types"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read LJ atom types.\n" ENDFB; } if(!amber7) { /* skip excluded atom counts */ p=skip_fortran(nAtom,12,p); /* skip NB param arrays */ p=skip_fortran(NTYPES*NTYPES,12,p); } /* read residue labels */ if(amber7) { p = findflag(buffer,"RESIDUE_LABEL","20a4"); } resn = Alloc(ResName,NRES); b=0; for(a=0;a<NRES;a++) { p=ncopy(cc,p,4); if(!sscanf(cc,"%s",resn[a])) resn[a][0]=0; if((++b)==20) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error reading residue labels"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read residue labels.\n" ENDFB; } /* read residue assignments */ if(amber7) { p = findflag(buffer,"RESIDUE_POINTER","10I8"); col=10; wid=8; } else { col=12; wid=6; } b=0; last_i=0; rc=0; for(a=0;a<NRES;a++) { p=ncopy(cc,p,wid); if(sscanf(cc,"%d",&at_i)) { if(last_i) for(aa=(last_i-1);aa<(at_i-1);aa++) { ai = atInfo+aa; strcpy(ai->resn,resn[a-1]); ai->resv=rc; sprintf(ai->resi,"%d",rc); } rc++; last_i=at_i; } if((++b)==col) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(last_i) for(aa=(last_i-1);aa<nAtom;aa++) { ai = atInfo+aa; strcpy(ai->resn,resn[NRES-1]); ai->resv=rc; sprintf(ai->resi,"%d",rc); } rc++; if(!ok) { ErrMessage("TOPStrToCoordSet","Error reading residues"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read residues.\n" ENDFB; } FreeP(resn); if(!amber7) { /* skip bond force constants */ p=skip_fortran(NUMBND,5,p); /* skip bond lengths */ p=skip_fortran(NUMBND,5,p); /* skip angle force constant */ p=skip_fortran(NUMANG,5,p); /* skip angle eq */ p=skip_fortran(NUMANG,5,p); /* skip dihedral force constant */ p=skip_fortran(NPTRA,5,p); /* skip dihedral periodicity */ p=skip_fortran(NPTRA,5,p); /* skip dihedral phases */ p=skip_fortran(NPTRA,5,p); /* skip SOLTYs */ p=skip_fortran(NATYP,5,p); /* skip LJ terms r12 */ p=skip_fortran((NTYPES*(NTYPES+1))/2,5,p); /* skip LJ terms r6 */ p=skip_fortran((NTYPES*(NTYPES+1))/2,5,p); } /* read bonds */ if(amber7) { p = findflag(buffer,"BONDS_INC_HYDROGEN","10I8"); col=10; wid=8; } else { col=12; wid=6; } nBond = NBONH + NBONA; bond=VLAlloc(BondType,nBond); bi = 0; b=0; c=0; i0=0; i1=0; for(a=0;a<3*NBONH;a++) { p=ncopy(cc,p,wid); i2=i1; i1=i0; if(!sscanf(cc,"%d",&i0)) ok=false; if((++c)==3) { c=0; bd=bond+bi; bd->index[0]=(abs(i2)/3); bd->index[1]=(abs(i1)/3); bd->order=1; bd->stereo=0; bd->id = bi+1; bi++; } if((++b)==col) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error hydrogen containing bonds"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read %d hydrogen containing bonds.\n",NBONH ENDFB; } if(amber7) { p = findflag(buffer,"BONDS_WITHOUT_HYDROGEN","10I8"); col=10; wid=8; } else { col=12; wid=6; } b=0; c=0; for(a=0;a<3*NBONA;a++) { p=ncopy(cc,p,wid); i2=i1; i1=i0; if(!sscanf(cc,"%d",&i0)) ok=false; if((++c)==3) { c=0; bd=bond+bi; bd->index[0]=(abs(i2)/3); bd->index[1]=(abs(i1)/3); bd->order=0; bd->stereo=0; bd->id = bi+1; bi++; } if((++b)==col) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error hydrogen free bonds"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read %d hydrogen free bonds.\n",NBONA ENDFB; } if(!amber7) { /* skip hydrogen angles */ p=skip_fortran(4*NTHETH,12,p); /* skip non-hydrogen angles */ p=skip_fortran(4*NTHETA,12,p); /* skip hydrogen dihedrals */ p=skip_fortran(5*NPHIH,12,p); /* skip non hydrogen dihedrals */ p=skip_fortran(5*NPHIA,12,p); /* skip nonbonded exclusions */ p=skip_fortran(NNB,12,p); /* skip hydrogen bonds ASOL */ p=skip_fortran(NPHB,5,p); /* skip hydrogen bonds BSOL */ p=skip_fortran(NPHB,5,p); /* skip HBCUT */ p=skip_fortran(NPHB,5,p); } /* read AMBER atom types */ if(amber7) { p = findflag(buffer,"AMBER_ATOM_TYPE","20a4"); } b=0; for(a=0;a<nAtom;a++) { p=ncopy(cc,p,4); ai=atInfo+a; if(!sscanf(cc,"%s",ai->textType)) ok=false; if((++b)==20) { b=0; p=nextline_top(p); } } if(b) p=nextline_top(p); if(!ok) { ErrMessage("TOPStrToCoordSet","Error reading atom types"); } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " TOPStr2CoordSet: read atom types.\n" ENDFB; } if(!amber7) { /* skip TREE classification */ p=skip_fortran(nAtom,20,p); /* skip tree joining information */ p=skip_fortran(nAtom,12,p); /* skip last atom rotated blah blah blah */ p=skip_fortran(nAtom,12,p); } if(IFBOX>0) { int IPTRES,NSPM,NSPSOL; if(amber7) { p = findflag(buffer,"SOLVENT_POINTERS","3I8"); wid=8; } else { wid=6; } p=ncopy(cc,p,wid); ok = ok && (sscanf(cc,"%d",&IPTRES)==1); p=ncopy(cc,p,wid); ok = ok && (sscanf(cc,"%d",&NSPM)==1); p=ncopy(cc,p,wid); ok = ok && (sscanf(cc,"%d",&NSPSOL)==1); p=nextline_top(p); if(amber7) { p = findflag(buffer,"ATOMS_PER_MOLECULE","10I8"); col=10; } else { col=12; } /* skip num atoms per box */ p=skip_fortran(NSPM,col,p); if(amber7) { p = findflag(buffer,"BOX_DIMENSIONS","5E16.8"); } wid=16; p=ncopy(cc,p,16); ok = ok && (sscanf(cc,"%f",&BETA)==1); p=ncopy(cc,p,16); ok = ok && (sscanf(cc,"%f",&BOX1)==1); p=ncopy(cc,p,16); ok = ok && (sscanf(cc,"%f",&BOX2)==1); p=ncopy(cc,p,16); ok = ok && (sscanf(cc,"%f",&BOX3)==1); if(ok) { if(!cset->PeriodicBox) cset->PeriodicBox=CrystalNew(); cset->PeriodicBox->Dim[0] = BOX1; cset->PeriodicBox->Dim[1] = BOX2; cset->PeriodicBox->Dim[2] = BOX3; if((BETA > 109.47) && (BETA < 109.48)) { cset->PeriodicBoxType=cCSet_Octahedral; cset->PeriodicBox->Angle[0]=(float)(2.0*acos(1.0/sqrt(3.0))*180.0/PI); cset->PeriodicBox->Angle[1]=(float)(2.0*acos(1.0/sqrt(3.0))*180.0/PI); cset->PeriodicBox->Angle[2]=(float)(2.0*acos(1.0/sqrt(3.0))*180.0/PI); } else if(BETA==60.0) { cset->PeriodicBox->Angle[0]=60.0; /* rhombic dodecahedron (from ptraj.c) */ cset->PeriodicBox->Angle[1]=90.0; cset->PeriodicBox->Angle[2]=60.0; } else { cset->PeriodicBox->Angle[0]=90.0; cset->PeriodicBox->Angle[1]=BETA; cset->PeriodicBox->Angle[2]=90.0; } CrystalUpdate(cset->PeriodicBox); /* CrystalDump(cset->PeriodicBox);*/ } /* skip periodic box */ p=nextline_top(p); } if(!amber7) { if(IFCAP>0) { p=nextline_top(p); p=nextline_top(p); p=nextline_top(p); } if(IFPERT>0) { /* skip perturbed bond atoms */ p=skip_fortran(2*NBPER,12,p); /* skip perturbed bond atom pointers */ p=skip_fortran(2*NBPER,12,p); /* skip perturbed angles */ p=skip_fortran(3*NGPER,12,p); /* skip perturbed angle pointers */ p=skip_fortran(2*NGPER,12,p); /* skip perturbed dihedrals */ p=skip_fortran(4*NDPER,12,p); /* skip perturbed dihedral pointers */ p=skip_fortran(2*NDPER,12,p); /* skip residue names */ p=skip_fortran(NRES,20,p); /* skip atom names */ p=skip_fortran(nAtom,20,p); /* skip atom symbols */ p=skip_fortran(nAtom,20,p); /* skip unused field */ p=skip_fortran(nAtom,5,p); /* skip perturbed flags */ p=skip_fortran(nAtom,12,p); /* skip LJ atom flags */ p=skip_fortran(nAtom,12,p); /* skip perturbed charges */ p=skip_fortran(nAtom,5,p); } if(IPOL>0) { /* skip atomic polarizabilities */ p=skip_fortran(nAtom,5,p); } if((IPOL>0) && (IFPERT>0)) { /* skip atomic polarizabilities */ p=skip_fortran(nAtom,5,p); } } /* for future reference %FLAG LES_NTYP %FORMAT(10I8) %FLAG LES_TYPE %FORMAT(10I8) %FLAG LES_FAC %FORMAT(5E16.8) %FLAG LES_CNUM %FORMAT(10I8) %FLAG LES_ID %FORMAT(10I8) Here is the additional information for LES topology formats: First, if NPARM ==1, LES entries are in topology (NPARM is the 10th entry in the initial list of control parameters); otherwise the standard format applies. So, with NPARM=1, you just need to read a few more things at the very end of topology file: LES_NTYP (format: I6) ... one number, number of LES types and four arrays: LES_TYPE (12I6) ... NATOM integer entries LES_FAC (E16.8) ... LES_NTYPxLES_NTYP float entries LES_CNUM (12I6) ... NATOM integer entries LES_ID (12I6) ... NATOM integer entries and that's it. Your parser must have skipped this information because it was at the end of the file. Maybe that's good enough. */ coord=VLAlloc(float,3*nAtom); f=coord; for(a=0;a<nAtom;a++) { *(f++)=0.0; *(f++)=0.0; *(f++)=0.0; ai = atInfo + a; ai->id = a+1; /* assign 1-based identifiers */ AtomInfoAssignParameters(ai); ai->color=AtomInfoGetColor(ai); for(c=0;c<cRepCnt;c++) { ai->visRep[c] = false; } ai->visRep[cRepLine] = auto_show_lines; /* show lines by default */ ai->visRep[cRepNonbonded] = auto_show_nonbonded; /* show lines by default */ } } if(ok) { cset->NIndex=nAtom; cset->Coord=coord; cset->TmpBond=bond; cset->NTmpBond=nBond; } else { if(cset) cset->fFree(cset); } if(atInfoPtr) *atInfoPtr = atInfo; return(cset); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadTOPStr(ObjectMolecule *I,char *TOPStr,int frame,int discrete) { CoordSet *cset = NULL; AtomInfoType *atInfo; int ok=true; int isNew = true; unsigned int nAtom = 0; if(!I) isNew=true; else isNew=false; if(ok) { if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { /* never */ atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset=ObjectMoleculeTOPStr2CoordSet(TOPStr,&atInfo); nAtom=cset->NIndex; } /* include coordinate set */ if(ok) { if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } else { ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_AllMask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; /* if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; */ if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); if(cset->Symmetry&&(!I->Symmetry)) { I->Symmetry=SymmetryCopy(cset->Symmetry); SymmetryAttemptGeneration(I->Symmetry,false,false); } if(I->CSTmpl) if(I->CSTmpl->fFree) I->CSTmpl->fFree(I->CSTmpl); I->CSTmpl = cset; /* save template coordinate set */ SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); } return(I); } ObjectMolecule *ObjectMoleculeLoadTOPFile(ObjectMolecule *obj,char *fname,int frame,int discrete) { ObjectMolecule *I=NULL; int ok=true; FILE *f; long size; char *buffer,*p; f=fopen(fname,"rb"); if(!f) ok=ErrMessage("ObjectMoleculeLoadTOPFile","Unable to open file!"); else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadTOPFile: Loading from %s.\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); I=ObjectMoleculeReadTOPStr(obj,buffer,frame,discrete); mfree(buffer); } return(I); } void ObjectMoleculeSculptClear(ObjectMolecule *I) { PRINTFD(FB_ObjectMolecule) " ObjectMoleculeSculptClear: entered.\n" ENDFD; if(I->Sculpt) SculptFree(I->Sculpt); I->Sculpt=NULL; } void ObjectMoleculeSculptImprint(ObjectMolecule *I,int state) { PRINTFD(FB_ObjectMolecule) " ObjectMoleculeUpdateSculpt: entered.\n" ENDFD; if(!I->Sculpt) I->Sculpt = SculptNew(); SculptMeasureObject(I->Sculpt,I,state); } float ObjectMoleculeSculptIterate(ObjectMolecule *I,int state,int n_cycle) { PRINTFD(FB_ObjectMolecule) " ObjectMoleculeIterateSculpt: entered.\n" ENDFD; if(I->Sculpt) { return SculptIterateObject(I->Sculpt,I,state,n_cycle); } else return 0.0F; } void ObjectMoleculeUpdateIDNumbers(ObjectMolecule *I) { int a; int max; AtomInfoType *ai; BondType *b; if(I->AtomCounter<0) { max=-1; ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(ai->id>max) max=ai->id; ai++; } I->AtomCounter=max+1; } ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(ai->id<0) ai->id=I->AtomCounter++; ai++; } if(I->BondCounter<0) { max=-1; b=I->Bond; for(a=0;a<I->NBond;a++) { if(b->id>max) max=b->id; b++; } I->BondCounter=max+1; } b=I->Bond; for(a=0;a<I->NBond;a++) { if(!b->id) b->id=I->BondCounter++; b++; } } CoordSet *ObjectMoleculePMO2CoordSet(CRaw *pmo,AtomInfoType **atInfoPtr,int *restart) { int nAtom,nBond; int a; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL,*ai; AtomInfoType068 *atInfo068 = NULL; AtomInfoType076 *atInfo076 = NULL; AtomInfoType083 *atInfo083 = NULL; BondType *bond=NULL; BondType068 *bond068=NULL; BondType083 *bond083=NULL; int ok=true; int auto_show_lines; int auto_show_nonbonded; int type,size; float *spheroid=NULL; float *spheroid_normal=NULL; int sph_info[2]; int version; auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); *restart=false; nAtom=0; nBond=0; if(atInfoPtr) atInfo = *atInfoPtr; type = RawGetNext(pmo,&size,&version); if(type!=cRaw_AtomInfo1) { ok=false; } else { /* read atoms */ PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loading atom info %d bytes = %8.3f\n",size,((float)size)/sizeof(AtomInfoType) ENDFD; if(version<66) { PRINTFB(FB_ObjectMolecule,FB_Errors) " ObjectMolecule: unsupported binary file (version %d). aborting.\n", version ENDFB; ok=false; } else if(version<69) { /* legacy atom format */ nAtom = size/sizeof(AtomInfoType068); atInfo068 = Alloc(AtomInfoType068,nAtom); ok = RawReadInto(pmo,cRaw_AtomInfo1,size,(char*)atInfo068); VLACheck(atInfo,AtomInfoType,nAtom); UtilExpandArrayElements(atInfo068,atInfo,nAtom, sizeof(AtomInfoType068),sizeof(AtomInfoType)); FreeP(atInfo068); } else if(version<77) { /* legacy atom format */ nAtom = size/sizeof(AtomInfoType076); atInfo076 = Alloc(AtomInfoType076,nAtom); ok = RawReadInto(pmo,cRaw_AtomInfo1,size,(char*)atInfo076); VLACheck(atInfo,AtomInfoType,nAtom); UtilExpandArrayElements(atInfo076,atInfo,nAtom, sizeof(AtomInfoType076),sizeof(AtomInfoType)); FreeP(atInfo076); } else if(version<84) { /* legacy atom format */ nAtom = size/sizeof(AtomInfoType083); atInfo083 = Alloc(AtomInfoType083,nAtom); ok = RawReadInto(pmo,cRaw_AtomInfo1,size,(char*)atInfo083); VLACheck(atInfo,AtomInfoType,nAtom); UtilExpandArrayElements(atInfo083,atInfo,nAtom, sizeof(AtomInfoType083),sizeof(AtomInfoType)); FreeP(atInfo083); } else { nAtom = size/sizeof(AtomInfoType); VLACheck(atInfo,AtomInfoType,nAtom); ok = RawReadInto(pmo,cRaw_AtomInfo1,size,(char*)atInfo); } } if(ok) { PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loading coordinates\n" ENDFD; coord = (float*)RawReadVLA(pmo,cRaw_Coords1,sizeof(float),5,false); if(!coord) ok=false; } type = RawGetNext(pmo,&size,&version); if(type==cRaw_SpheroidInfo1) { PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loading spheroid\n" ENDFD; ok = RawReadInto(pmo,cRaw_SpheroidInfo1,sizeof(int)*2,(char*)sph_info); if(ok) { PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loading spheroid size %d nsph %d\n",sph_info[0],sph_info[1] ENDFD; spheroid = (float*)RawReadPtr(pmo,cRaw_Spheroid1,&size); if(!spheroid) ok=false; PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loaded spheroid %p size %d \n", (void*)spheroid,size ENDFD; } if(ok) { spheroid_normal = (float*)RawReadPtr(pmo,cRaw_SpheroidNormals1,&size); if(!spheroid_normal) ok=false; } PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loaded spheroid %p size %d \n", (void*)spheroid_normal,size ENDFD; } if(ok) type = RawGetNext(pmo,&size,&version); if(ok) { PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: loading bonds\n" ENDFD; if(type!=cRaw_Bonds1) { ok=false; } else { if(ok) { /* legacy bond format */ if(version<66) { PRINTFB(FB_ObjectMolecule,FB_Errors) " ObjectMolecule: unsupported binary file (version %d). aborting.\n", version ENDFB; ok=false; } else if(version<69) { /* legacy atom format */ nBond = size/sizeof(BondType068); bond068 = Alloc(BondType068,nBond); ok = RawReadInto(pmo,cRaw_Bonds1,nBond*sizeof(BondType068),(char*)bond068); bond=VLAlloc(BondType,nBond); UtilExpandArrayElements(bond068,bond,nBond, sizeof(BondType068),sizeof(BondType)); FreeP(bond068); for(a=0;a<nBond;a++) bond[a].id=-1; /* initialize identifiers */ } else if(version<84) { nBond = size/sizeof(BondType083); bond083 = Alloc(BondType083,nBond); ok = RawReadInto(pmo,cRaw_Bonds1,nBond*sizeof(BondType083),(char*)bond083); bond=VLAlloc(BondType,nBond); UtilExpandArrayElements(bond083,bond,nBond, sizeof(BondType083),sizeof(BondType)); FreeP(bond083); } else { bond=(BondType*)RawReadVLA(pmo,cRaw_Bonds1,sizeof(BondType),5,false); nBond = VLAGetSize(bond); } PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoordSet: found %d bonds\n",nBond ENDFD; if(Feedback(FB_ObjectMolecule,FB_Debugging)) { for(a=0;a<nBond;a++) printf(" ObjectMoleculeConnect: bond %d ind0 %d ind1 %d order %d\n", a,bond[a].index[0],bond[a].index[1],bond[a].order); } } } } if(ok) { ai=atInfo; for(a=0;a<nAtom;a++) { ai->selEntry=0; ai++; } cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->NTmpBond=nBond; cset->TmpBond=bond; if(spheroid) { cset->Spheroid=spheroid; cset->SpheroidNormal=spheroid_normal; cset->SpheroidSphereSize=sph_info[0]; cset->NSpheroid = sph_info[1]; } } else { VLAFreeP(bond); VLAFreeP(coord); FreeP(spheroid); FreeP(spheroid_normal); } if(atInfoPtr) *atInfoPtr = atInfo; if(ok) { type = RawGetNext(pmo,&size,&version); if(type==cRaw_AtomInfo1) *restart=true; } return(cset); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadPMO(ObjectMolecule *I,CRaw *pmo,int frame,int discrete) { CoordSet *cset = NULL; AtomInfoType *atInfo; int ok=true; int isNew = true; unsigned int nAtom = 0; int restart =false; int repeatFlag = true; int successCnt = 0; while(repeatFlag) { repeatFlag = false; if(!I) isNew=true; else isNew=false; if(ok) { if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset = ObjectMoleculePMO2CoordSet(pmo,&atInfo,&restart); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } if(cset) nAtom=cset->NIndex; else ok=false; } /* include coordinate set */ if(ok) { if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(!isNew) { ObjectMoleculeMerge(I,atInfo,cset,true,cAIC_AllMask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); if(cset->Symmetry&&(!I->Symmetry)) { I->Symmetry=SymmetryCopy(cset->Symmetry); SymmetryAttemptGeneration(I->Symmetry,false,false); } SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); successCnt++; if(successCnt>1) { if(successCnt==2){ PRINTFB(FB_ObjectMolecule,FB_Actions) " ObjectMolReadPMO: read model %d\n",1 ENDFB; } PRINTFB(FB_ObjectMolecule,FB_Actions) " ObjectMolReadPMO: read model %d\n",successCnt ENDFB; } } if(restart) { repeatFlag=true; frame=frame+1; restart=false; } } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadPMOFile(ObjectMolecule *obj,char *fname,int frame,int discrete) { ObjectMolecule *I=NULL; int ok=true; CRaw *raw; raw = RawOpenRead(fname); if(!raw) ok=ErrMessage("ObjectMoleculeLoadPMOFile","Unable to open file!"); else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadPMOFile: Loading from %s.\n",fname ENDFB; I=ObjectMoleculeReadPMO(obj,raw,frame,discrete); RawFree(raw); } return(I); } /*========================================================================*/ int ObjectMoleculeMultiSave(ObjectMolecule *I,char *fname,int state,int append) { /* version 1 writes atominfo, coords, spheroid, bonds */ CRaw *raw = NULL; int ok=true; int a,c,a1,a2,b1,b2; BondType *b; CoordSet *cs; BondType *bondVLA = NULL; AtomInfoType *aiVLA = NULL; int start,stop; int nBond; int sph_info[2]; PRINTFD(FB_ObjectMolecule) " ObjectMoleculeMultiSave-Debug: entered \"%s\" state=%d\n",fname,state ENDFD; if(append) { raw = RawOpenWrite(fname); } else { raw = RawOpenAppend(fname); } if(raw) { aiVLA = VLAMalloc(1000,sizeof(AtomInfoType),5,true); bondVLA = VLAlloc(BondType,4000); if(state<0) { start=0; stop=I->NCSet; } else { start=state; if(start<0) start=0; stop=state+1; if(stop>I->NCSet) stop=I->NCSet; } for(a=start;a<stop;a++) { PRINTFD(FB_ObjectMolecule) " ObjectMMSave-Debug: state %d\n",a ENDFD; cs=I->CSet[a]; if(cs) { VLACheck(aiVLA,AtomInfoType,cs->NIndex); nBond=0; /* write atoms */ for(c=0;c<cs->NIndex;c++) { a1 = cs->IdxToAtm[c]; /* always valid */ aiVLA[c]=I->AtomInfo[a1]; } if(ok) ok = RawWrite(raw,cRaw_AtomInfo1,sizeof(AtomInfoType)*cs->NIndex,0,(char*)aiVLA); /* write coords */ if(ok) ok = RawWrite(raw,cRaw_Coords1,sizeof(float)*3*cs->NIndex,0,(char*)cs->Coord); /* write spheroid (if one exists) */ if(cs->Spheroid&&cs->SpheroidNormal) { sph_info[0]=cs->SpheroidSphereSize; sph_info[1]=cs->NSpheroid; if(ok) ok = RawWrite(raw,cRaw_SpheroidInfo1,sizeof(int)*2,0,(char*)sph_info); if(ok) ok = RawWrite(raw,cRaw_Spheroid1,sizeof(float)*cs->NSpheroid,0,(char*)cs->Spheroid); if(ok) ok = RawWrite(raw,cRaw_SpheroidNormals1,sizeof(float)*3*cs->NSpheroid,0,(char*)cs->SpheroidNormal); PRINTFD(FB_ObjectMolecule) " ObjectMolPMO2CoorSet: saved spheroid size %d %d\n",cs->SpheroidSphereSize,cs->NSpheroid ENDFD; } /* write bonds */ b=I->Bond; for(c=0;c<I->NBond;c++) { b1 = b->index[0]; b2 = b->index[1]; if(I->DiscreteFlag) { if((cs==I->DiscreteCSet[b1])&&(cs==I->DiscreteCSet[b2])) { a1=I->DiscreteAtmToIdx[b1]; a2=I->DiscreteAtmToIdx[b2]; } else { a1=-1; a2=-1; } } else { a1=cs->AtmToIdx[b1]; a2=cs->AtmToIdx[b2]; } if((a1>=0)&&(a2>=0)) { nBond++; VLACheck(bondVLA,BondType,nBond); bondVLA[nBond-1]=*b; bondVLA[nBond-1].index[0] = a1; bondVLA[nBond-1].index[1] = a2; } b++; } if(ok) ok = RawWrite(raw,cRaw_Bonds1,sizeof(BondType)*nBond,0,(char*)bondVLA); } } } if(raw) RawFree(raw); VLAFreeP(aiVLA); VLAFreeP(bondVLA); return(ok); } /*========================================================================*/ int ObjectMoleculeGetPhiPsi(ObjectMolecule *I,int ca,float *phi,float *psi,int state) { int np=-1; int cm=-1; int c=-1; int n=-1; int result = false; AtomInfoType *ai; int n0,at; float v_ca[3]; float v_n[3]; float v_c[3]; float v_cm[3]; float v_np[3]; ai=I->AtomInfo; if((ai[ca].name[0]=='C')&&(ai[ca].name[1]=='A')) { ObjectMoleculeUpdateNeighbors(I); /* find C */ n0 = I->Neighbor[ca]+1; while(I->Neighbor[n0]>=0) { at = I->Neighbor[n0]; if((ai[at].name[0]=='C')&&(ai[at].name[1]==0)) { c=at; break; } n0+=2; } /* find N */ n0 = I->Neighbor[ca]+1; while(I->Neighbor[n0]>=0) { at = I->Neighbor[n0]; if((ai[at].name[0]=='N')&&(ai[at].name[1]==0)) { n=at; break; } n0+=2; } /* find NP */ if(c>=0) { n0 = I->Neighbor[c]+1; while(I->Neighbor[n0]>=0) { at = I->Neighbor[n0]; if((ai[at].name[0]=='N')&&(ai[at].name[1]==0)) { np=at; break; } n0+=2; } } /* find CM */ if(n>=0) { n0 = I->Neighbor[n]+1; while(I->Neighbor[n0]>=0) { at = I->Neighbor[n0]; if((ai[at].name[0]=='C')&&(ai[at].name[1]==0)) { cm=at; break; } n0+=2; } } if((ca>=0)&&(np>=0)&&(c>=0)&&(n>=0)&&(cm>=0)) { if(ObjectMoleculeGetAtomVertex(I,state,ca,v_ca)&& ObjectMoleculeGetAtomVertex(I,state,n,v_n)&& ObjectMoleculeGetAtomVertex(I,state,c,v_c)&& ObjectMoleculeGetAtomVertex(I,state,cm,v_cm)&& ObjectMoleculeGetAtomVertex(I,state,np,v_np)) { (*phi)=rad_to_deg(get_dihedral3f(v_c,v_ca,v_n,v_cm)); (*psi)=rad_to_deg(get_dihedral3f(v_np,v_c,v_ca,v_n)); result=true; } } } return(result); } /*========================================================================*/ int ObjectMoleculeCheckBondSep(ObjectMolecule *I,int a0,int a1,int dist) { int result = false; int n0; int stack[MAX_BOND_DIST+1]; int history[MAX_BOND_DIST+1]; int depth=0; int distinct; int a; if(dist>MAX_BOND_DIST) return false; ObjectMoleculeUpdateNeighbors(I); PRINTFD(FB_ObjectMolecule) " CBS-Debug: %s %d %d %d\n",I->Obj.Name,a0,a1,dist ENDFD; depth = 1; history[depth]=a0; stack[depth] = I->Neighbor[a0]+1; /* go to first neighbor */ while(depth) { /* keep going until we've traversed tree */ while(I->Neighbor[stack[depth]]>=0) /* end of branches? go back up one bond */ { n0 = I->Neighbor[stack[depth]]; /* get current neighbor index */ stack[depth]+=2; /* set up next neighbor */ distinct=true; /* check to see if current candidate is distinct from ancestors */ for(a=1;a<depth;a++) { if(history[a]==n0) distinct=false; } if(distinct) { if(depth<dist) { /* are not yet at the proper distance? */ if(distinct) { depth++; stack[depth] = I->Neighbor[n0]+1; /* then keep moving outward */ history[depth] = n0; } } else if(n0==a1) /* otherwise, see if we have a match */ result = true; } } depth--; } PRINTFD(FB_ObjectMolecule) " CBS-Debug: result %d\n",result ENDFD; return result; } /*========================================================================*/ void ObjectGotoState(ObjectMolecule *I,int state) { if((I->NCSet>1)||(!SettingGet(cSetting_static_singletons))) { if(state>I->NCSet) state = I->NCSet-1; if(state<0) state = I->NCSet-1; SceneSetFrame(0,state); } } /*========================================================================*/ CSetting **ObjectMoleculeGetSettingHandle(ObjectMolecule *I,int state) { if(state<0) { return(&I->Obj.Setting); } else if(state<I->NCSet) { if(I->CSet[state]) { return(&I->CSet[state]->Setting); } else { return(NULL); } } else { return(NULL); } } /*========================================================================*/ int ObjectMoleculeSetStateTitle(ObjectMolecule *I,int state,char *text) { int result=false; if(state<0) state=I->NCSet-1; if(state>=I->NCSet) { PRINTFB(FB_ObjectMolecule,FB_Errors) "Error: invalid state %d\n",state +1 ENDFB; } else if(!I->CSet[state]) { PRINTFB(FB_ObjectMolecule,FB_Errors) "Error: empty state %d\n",state +1 ENDFB; } else { UtilNCopy(I->CSet[state]->Name,text,sizeof(WordType)); result=true; } return(result); } /*========================================================================*/ char *ObjectMoleculeGetStateTitle(ObjectMolecule *I,int state) { char *result=NULL; if(state<0) state=I->NCSet-1; if(state>=I->NCSet) { PRINTFB(FB_ObjectMolecule,FB_Errors) "Error: invalid state %d\n",state +1 ENDFB; } else if(!I->CSet[state]) { PRINTFB(FB_ObjectMolecule,FB_Errors) "Error: empty state %d\n",state +1 ENDFB; } else { result = I->CSet[state]->Name; } return(result); } /*========================================================================*/ void ObjectMoleculeRenderSele(ObjectMolecule *I,int curState,int sele) { CoordSet *cs; int a,at; if(PMGUI) { if(curState>=0) { if(curState<I->NCSet) { if(I->CSet[curState]) { cs=I->CSet[curState]; for(a=0;a<cs->NIndex;a++) { at=cs->IdxToAtm[a]; /* should work for both discrete and non-discrete objects */ if(SelectorIsMember(I->AtomInfo[at].selEntry,sele)) glVertex3fv(cs->Coord+3*a); } } } else if(SettingGet(cSetting_static_singletons)) { if(I->NCSet==1) { cs=I->CSet[0]; if(cs) { for(a=0;a<cs->NIndex;a++) { at=cs->IdxToAtm[a]; /* should work for both discrete and non-discrete objects */ if(SelectorIsMember(I->AtomInfo[at].selEntry,sele)) glVertex3fv(cs->Coord+3*a); } } } } } else { /* all states */ for(curState=0;curState<I->NCSet;curState++) { if(I->CSet[curState]) { cs=I->CSet[curState]; for(a=0;a<cs->NIndex;a++) { at=cs->IdxToAtm[a]; /* should work for both discrete and non-discrete objects */ if(SelectorIsMember(I->AtomInfo[at].selEntry,sele)) glVertex3fv(cs->Coord+3*a); } } } } } } /*========================================================================*/ CoordSet *ObjectMoleculeXYZStr2CoordSet(char *buffer,AtomInfoType **atInfoPtr) { char *p; int nAtom; int a,c; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL,*ai; char cc[MAXLINELEN]; int atomCount; BondType *bond=NULL; int nBond=0; int b1,b2; WordType tmp_name; int auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); int auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); BondType *ii; p=buffer; nAtom=0; atInfo = *atInfoPtr; p=ncopy(cc,p,6); if(!sscanf(cc,"%d",&nAtom)) nAtom=0; p=nskip(p,2); p=ncopy(tmp_name,p,sizeof(WordType)-1); p=nextline_top(p); coord=VLAlloc(float,3*nAtom); if(atInfo) VLACheck(atInfo,AtomInfoType,nAtom); nBond=0; bond=VLAlloc(BondType,6*nAtom); ii=bond; PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeReadXYZ: Found %i atoms...\n",nAtom ENDFB; a=0; atomCount=0; while(*p) { ai=atInfo+atomCount; p=ncopy(cc,p,6); if(!sscanf(cc,"%d",&ai->id)) break; p=nskip(p,2);/* to 12 */ p=ncopy(cc,p,3); if(!sscanf(cc,"%s",ai->name)) ai->name[0]=0; ai->alt[0]=0; strcpy(ai->resn,"UNK"); ai->chain[0] = 0; ai->resv=atomCount+1; sprintf(ai->resi,"%d",ai->resv); p=ncopy(cc,p,12); sscanf(cc,"%f",coord+a); p=ncopy(cc,p,12); sscanf(cc,"%f",coord+(a+1)); p=ncopy(cc,p,12); sscanf(cc,"%f",coord+(a+2)); ai->q=1.0; ai->b=0.0; ai->segi[0]=0; ai->elem[0]=0; /* let atom info guess/infer atom type */ for(c=0;c<cRepCnt;c++) { ai->visRep[c] = false; } ai->visRep[cRepLine] = auto_show_lines; /* show lines by default */ ai->visRep[cRepNonbonded] = auto_show_nonbonded; /* show lines by default */ p=ncopy(cc,p,6); sscanf(cc,"%d",&ai->customType); /* in the absense of external tinker information, assume hetatm */ ai->hetatm=1; AtomInfoAssignParameters(ai); ai->color=AtomInfoGetColor(ai); b1 = atomCount; for(c=0;c<6;c++) { p=ncopy(cc,p,6); if (!cc[0]) break; if(!sscanf(cc,"%d",&b2)) break; if(b1<(b2-1)) { nBond++; ii->index[0] = b1; ii->index[1] = b2-1; ii->order = 1; /* missing bond order information */ ii->stereo = 0; ii->id = -1; /* no serial number */ } } PRINTFD(FB_ObjectMolecule) " ObjectMolecule-DEBUG: %s %s %s %s %8.3f %8.3f %8.3f %6.2f %6.2f %s\n", ai->name,ai->resn,ai->resi,ai->chain, *(coord+a),*(coord+a+1),*(coord+a+2),ai->b,ai->q, ai->segi ENDFD; a+=3; atomCount++; if(atomCount>=nAtom) break; p=nextline_top(p); } PRINTFB(FB_ObjectMolecule,FB_Blather) " XYZStr2CoordSet: Read %d bonds.\n",nBond ENDFB; cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->TmpBond=bond; cset->NTmpBond=nBond; strcpy(cset->Name,tmp_name); if(atInfoPtr) *atInfoPtr = atInfo; return(cset); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadXYZStr(ObjectMolecule *I,char *PDBStr,int frame,int discrete) { CoordSet *cset = NULL; AtomInfoType *atInfo; int ok=true; int isNew = true; unsigned int nAtom = 0; if(!I) isNew=true; else isNew=false; if(ok) { if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset=ObjectMoleculeXYZStr2CoordSet(PDBStr,&atInfo); nAtom=cset->NIndex; } /* include coordinate set */ if(ok) { if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } else { ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_IDMask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); if(cset->Symmetry&&(!I->Symmetry)) { I->Symmetry=SymmetryCopy(cset->Symmetry); SymmetryAttemptGeneration(I->Symmetry,false,false); } SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadXYZFile(ObjectMolecule *obj,char *fname,int frame,int discrete) { ObjectMolecule *I=NULL; int ok=true; FILE *f; long size; char *buffer,*p; f=fopen(fname,"rb"); if(!f) ok=ErrMessage("ObjectMoleculeLoadXYZFile","Unable to open file!"); else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadXYZFile: Loading from %s.\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); I=ObjectMoleculeReadXYZStr(obj,buffer,frame,discrete); mfree(buffer); } return(I); } /*========================================================================*/ int ObjectMoleculeAreAtomsBonded(ObjectMolecule *I,int i0,int i1) { int result=false; int a; BondType *b; b=I->Bond; for (a=0;a<I->NBond;a++) { if(i0==b->index[0]) { if(i1==b->index[1]) { result=true; break; } } if(i1==b->index[0]) { if(i0==b->index[1]) { result=true; break; } } b++; } return(result); } /*========================================================================*/ void ObjectMoleculeRenameAtoms(ObjectMolecule *I,int force) { AtomInfoType *ai; int a; if(force) { ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) (ai++)->name[0]=0; } AtomInfoUniquefyNames(NULL,0,I->AtomInfo,I->NAtom); } /*========================================================================*/ void ObjectMoleculeAddSeleHydrogens(ObjectMolecule *I,int sele) { int a,b; int n,nn; CoordSet *cs; CoordSet *tcs; int seleFlag=false; AtomInfoType *ai,*nai,fakeH; int repeatFlag=false; int nH; int *index; float v[3],v0[3]; float d; UtilZeroMem(&fakeH,sizeof(AtomInfoType)); fakeH.protons=1; ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(SelectorIsMember(ai->selEntry,sele)) { seleFlag=true; break; } ai++; } if(seleFlag) { if(!ObjectMoleculeVerifyChemistry(I)) { ErrMessage(" AddHydrogens","missing chemical geometry information."); } else if(I->DiscreteFlag) { ErrMessage(" AddHydrogens","can't modify a discrete object."); } else { repeatFlag=true; while(repeatFlag) { repeatFlag=false; nH = 0; ObjectMoleculeUpdateNeighbors(I); nai = (AtomInfoType*)VLAMalloc(1000,sizeof(AtomInfoType),1,true); ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(SelectorIsMember(ai->selEntry,sele)) { n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(nn<ai->valence) { VLACheck(nai,AtomInfoType,nH); UtilNCopy((nai+nH)->elem,"H",2); (nai+nH)->geom=cAtomInfoSingle; (nai+nH)->valence=1; (nai+nH)->temp1 = a; /* borrowing this field temporarily */ ObjectMoleculePrepareAtom(I,a,nai+nH); nH++; } } ai++; } if(nH) { repeatFlag=true; cs = CoordSetNew(); cs->Coord = VLAlloc(float,nH*3); cs->NIndex=nH; index = Alloc(int,nH); for(a=0;a<nH;a++) { index[a] = (nai+a)->temp1; } if(cs->fEnumIndices) cs->fEnumIndices(cs); cs->TmpLinkBond = VLAlloc(BondType,nH); for(a=0;a<nH;a++) { cs->TmpLinkBond[a].index[0] = (nai+a)->temp1; cs->TmpLinkBond[a].index[1] = a; cs->TmpLinkBond[a].order = 1; cs->TmpLinkBond[a].stereo = 0; cs->TmpLinkBond[a].id = -1; } cs->NTmpLinkBond = nH; AtomInfoUniquefyNames(I->AtomInfo,I->NAtom,nai,nH); ObjectMoleculeMerge(I,nai,cs,false,cAIC_AllMask); /* will free nai and cs->TmpLinkBond */ ObjectMoleculeExtendIndices(I); ObjectMoleculeUpdateNeighbors(I); for(b=0;b<I->NCSet;b++) { /* add coordinate into the coordinate set */ tcs = I->CSet[b]; if(tcs) { for(a=0;a<nH;a++) { ObjectMoleculeGetAtomVertex(I,b,index[a],v0); ObjectMoleculeFindOpenValenceVector(I,b,index[a],v,NULL); d = AtomInfoGetBondLength(I->AtomInfo+index[a],&fakeH); scale3f(v,d,v); add3f(v0,v,cs->Coord+3*a); } CoordSetMerge(tcs,cs); } } FreeP(index); if(cs->fFree) cs->fFree(cs); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); } else VLAFreeP(nai); } } } } /*========================================================================*/ void ObjectMoleculeFuse(ObjectMolecule *I,int index0,ObjectMolecule *src,int index1,int mode) { int a,b; AtomInfoType *ai0,*ai1,*nai; int n,nn; int at0=-1; int at1=-1; int a0,a1; int hydr1=-1; int anch1=-1; int ca0,ch0; BondType *b0,*b1; float *backup = NULL; float d,*f0,*f1; float va0[3],vh0[3],va1[3],vh1[3]; float x0[3],y0[3],z0[3]; float x1[3],y1[3],z1[3]; float x[3],y[3],z[3]; float t[3],t2[3]; CoordSet *cs=NULL,*scs=NULL; int state1 = 0; CoordSet *tcs; int edit=1; OrthoLineType sele1,sele2,s1,s2; ObjectMoleculeUpdateNeighbors(I); ObjectMoleculeUpdateNeighbors(src); /* make sure each link point has only one neighbor */ ai0=I->AtomInfo; ai1=src->AtomInfo; switch(mode) { case 0: /* fusing by replacing hydrogens */ n = I->Neighbor[index0]; nn = I->Neighbor[n++]; if(nn==1) at0 = I->Neighbor[n]; n = src->Neighbor[index1]; nn = src->Neighbor[n++]; if(nn==1) at1 = src->Neighbor[n]; if(src->NCSet) { scs = src->CSet[state1]; anch1 = scs->AtmToIdx[at1]; hydr1 = scs->AtmToIdx[index1]; } break; case 1: /* fuse merely by drawing a bond */ at0 = index0; at1 = index1; if(src->NCSet) { scs = src->CSet[state1]; anch1 = scs->AtmToIdx[at1]; } break; } if((at0>=0)&&(at1>=0)&&scs&&(anch1>=0)) { /* have anchors and source coordinate set */ nai=(AtomInfoType*)VLAMalloc(src->NAtom,sizeof(AtomInfoType),1,true); /* copy atoms and atom info into a 1:1 direct mapping */ cs = CoordSetNew(); cs->Coord = VLAlloc(float,scs->NIndex*3); cs->NIndex = scs->NIndex; for(a=0;a<scs->NIndex;a++) { copy3f(scs->Coord+a*3,cs->Coord+a*3); a1 = scs->IdxToAtm[a]; *(nai+a) = *(ai1+a1); (nai+a)->selEntry=0; /* avoid duplicating selection references -> leads to hangs ! */ (nai+a)->temp1=0; /* clear marks */ } nai[at1].temp1=2; /* mark the connection point */ /* copy internal bond information*/ cs->TmpBond = VLAlloc(BondType,src->NBond); b1 = src->Bond; b0 = cs->TmpBond; cs->NTmpBond=0; for(a=0;a<src->NBond;a++) { a0 = scs->AtmToIdx[b1->index[0]]; a1 = scs->AtmToIdx[b1->index[1]]; if((a0>=0)&&(a1>=0)) { *b0=*b1; b0->index[0] = a0; b0->index[1] = a1; b0++; cs->NTmpBond++; } b1++; } backup = Alloc(float,cs->NIndex*3); /* make untransformed copy of coordinate set */ for(a=0;a<cs->NIndex;a++) { copy3f(cs->Coord+a*3,backup+a*3); } switch(mode) { case 0: nai[hydr1].deleteFlag=true; I->AtomInfo[index0].deleteFlag=true; copy3f(backup+3*anch1,va1); copy3f(backup+3*hydr1,vh1); subtract3f(va1,vh1,x1); /* note reverse dir from above */ get_system1f3f(x1,y1,z1); break; case 1: copy3f(backup+3*anch1,va1); ObjectMoleculeFindOpenValenceVector(src,state1,at1,x1,NULL); scale3f(x1,-1.0F,x1); get_system1f3f(x1,y1,z1); break; } /* set up the linking bond */ cs->TmpLinkBond = VLAlloc(BondType,1); cs->NTmpLinkBond = 1; cs->TmpLinkBond->index[0] = at0; cs->TmpLinkBond->index[1] = anch1; cs->TmpLinkBond->order = 1; cs->TmpLinkBond->stereo = 0; cs->TmpLinkBond->id = -1; if(cs->fEnumIndices) cs->fEnumIndices(cs); d = AtomInfoGetBondLength(ai0+at0,ai1+at1); AtomInfoUniquefyNames(I->AtomInfo,I->NAtom,nai,cs->NIndex); /* set up tags which will enable use to continue editing bond */ if(edit) { for(a=0;a<I->NAtom;a++) { I->AtomInfo[a].temp1=0; } I->AtomInfo[at0].temp1=1; } ObjectMoleculeMerge(I,nai,cs,false,cAIC_AllMask); /* will free nai, cs->TmpBond and cs->TmpLinkBond */ ObjectMoleculeExtendIndices(I); ObjectMoleculeUpdateNeighbors(I); for(a=0;a<I->NCSet;a++) { /* add coordinate into the coordinate set */ tcs = I->CSet[a]; if(tcs) { switch(mode) { case 0: ca0 = tcs->AtmToIdx[at0]; /* anchor */ ch0 = tcs->AtmToIdx[index0]; /* hydrogen */ if((ca0>=0)&&(ch0>=0)) { copy3f(tcs->Coord+3*ca0,va0); copy3f(tcs->Coord+3*ch0,vh0); subtract3f(vh0,va0,x0); get_system1f3f(x0,y0,z0); } break; case 1: ca0 = tcs->AtmToIdx[at0]; /* anchor */ if(ca0>=0) { ObjectMoleculeFindOpenValenceVector(I,a,at0,x0,NULL); copy3f(tcs->Coord+3*ca0,va0); get_system1f3f(x0,y0,z0); } break; } scale3f(x0,d,t2); add3f(va0,t2,t2); f0=backup; f1=cs->Coord; for(b=0;b<cs->NIndex;b++) { /* brute force transformation */ subtract3f(f0,va1,t); scale3f(x0,dot_product3f(t,x1),x); scale3f(y0,dot_product3f(t,y1),y); scale3f(z0,dot_product3f(t,z1),z); add3f(x,y,y); add3f(y,z,f1); add3f(t2,f1,f1); f0+=3; f1+=3; } CoordSetMerge(tcs,cs); } } switch(mode) { case 0: ObjectMoleculePurge(I); break; } ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); if(edit) { /* edit the resulting bond */ at0=-1; at1=-1; for(a=0;a<I->NAtom;a++) { if(I->AtomInfo[a].temp1==1) at0=a; if(I->AtomInfo[a].temp1==2) at1=a; } if((at0>=0)&&(at1>=0)) { sprintf(sele1,"%s`%d",I->Obj.Name,at1+1); /* points outward... */ sprintf(sele2,"%s`%d",I->Obj.Name,at0+1); SelectorGetTmp(sele1,s1); SelectorGetTmp(sele2,s2); EditorSelect(s1,s2,NULL,NULL,false,true,true); SelectorFreeTmp(s1); SelectorFreeTmp(s2); } } } if(cs) if(cs->fFree) cs->fFree(cs); FreeP(backup); } /*========================================================================*/ int ObjectMoleculeVerifyChemistry(ObjectMolecule *I) { int result=false; AtomInfoType *ai; int a; int flag; ai=I->AtomInfo; flag=true; for(a=0;a<I->NAtom;a++) { if(!ai->chemFlag) { flag=false; } ai++; } if(!flag) { if(I->CSet[0]) { /* right now this stuff is locked to state 0 */ ObjectMoleculeInferChemFromBonds(I,0); ObjectMoleculeInferChemFromNeighGeom(I,0); ObjectMoleculeInferHBondFromChem(I); /* ObjectMoleculeInferChemForProtein(I,0);*/ } flag=true; ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(!ai->chemFlag) { flag=false; break; } ai++; } } if(flag) result=true; return(result); } /*========================================================================*/ void ObjectMoleculeAttach(ObjectMolecule *I,int index,AtomInfoType *nai) { int a; AtomInfoType *ai; int n,nn; float v[3],v0[3],d; CoordSet *cs; ObjectMoleculeUpdateNeighbors(I); ai=I->AtomInfo+index; n = I->Neighbor[index]; nn = I->Neighbor[n++]; cs = CoordSetNew(); cs->Coord = VLAlloc(float,3); cs->NIndex=1; cs->TmpLinkBond = VLAlloc(BondType,1); cs->NTmpLinkBond = 1; cs->TmpLinkBond->index[0]=index; cs->TmpLinkBond->index[1]=0; cs->TmpLinkBond->order=1; cs->TmpLinkBond->stereo=0; cs->TmpLinkBond->id = -1; if(cs->fEnumIndices) cs->fEnumIndices(cs); ObjectMoleculePrepareAtom(I,index,nai); d = AtomInfoGetBondLength(ai,nai); ObjectMoleculeMerge(I,nai,cs,false,cAIC_AllMask); /* will free nai and cs->TmpLinkBond */ ObjectMoleculeExtendIndices(I); ObjectMoleculeUpdateNeighbors(I); for(a=0;a<I->NCSet;a++) { /* add atom to each coordinate set */ if(I->CSet[a]) { ObjectMoleculeGetAtomVertex(I,a,index,v0); ObjectMoleculeFindOpenValenceVector(I,a,index,v,NULL); scale3f(v,d,v); add3f(v0,v,cs->Coord); CoordSetMerge(I->CSet[a],cs); } } ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); if(cs->fFree) cs->fFree(cs); } /*========================================================================*/ int ObjectMoleculeFillOpenValences(ObjectMolecule *I,int index) { int a; AtomInfoType *ai,*nai; int n,nn; int result=0; int flag = true; float v[3],v0[3],d; CoordSet *cs; if((index>=0)&&(index<=I->NAtom)) { while(1) { ObjectMoleculeUpdateNeighbors(I); ai=I->AtomInfo+index; n = I->Neighbor[index]; nn = I->Neighbor[n++]; if((nn>=ai->valence)||(!flag)) break; flag=false; cs = CoordSetNew(); cs->Coord = VLAlloc(float,3); cs->NIndex=1; cs->TmpLinkBond = VLAlloc(BondType,1); cs->NTmpLinkBond = 1; cs->TmpLinkBond->index[0]=index; cs->TmpLinkBond->index[1]=0; cs->TmpLinkBond->order=1; cs->TmpLinkBond->stereo=0; cs->TmpLinkBond->id = -1; if(cs->fEnumIndices) cs->fEnumIndices(cs); nai = (AtomInfoType*)VLAMalloc(1,sizeof(AtomInfoType),1,true); UtilNCopy(nai->elem,"H",2); nai->geom=cAtomInfoSingle; nai->valence=1; ObjectMoleculePrepareAtom(I,index,nai); d = AtomInfoGetBondLength(ai,nai); ObjectMoleculeMerge(I,nai,cs,false,cAIC_AllMask); /* will free nai and cs->TmpLinkBond */ ObjectMoleculeExtendIndices(I); ObjectMoleculeUpdateNeighbors(I); for(a=0;a<I->NCSet;a++) { /* add atom to each coordinate set */ if(I->CSet[a]) { ObjectMoleculeGetAtomVertex(I,a,index,v0); ObjectMoleculeFindOpenValenceVector(I,a,index,v,NULL); scale3f(v,d,v); add3f(v0,v,cs->Coord); CoordSetMerge(I->CSet[a],cs); } } if(cs->fFree) cs->fFree(cs); result++; flag=true; } } ObjectMoleculeUpdateIDNumbers(I); return(result); } #define MaxOcc 100 /*========================================================================*/ static int get_planer_normal(ObjectMolecule *I,int state,int index,float *normal) { /* NOTE assumes neighbors are defined */ int found = false; int nOcc = 0; float occ[MaxOcc*3]; AtomInfoType *ai = I->AtomInfo + index; int n,a1; float v0[3],v1[3],v2[3],n0[3]; if(ObjectMoleculeGetAtomVertex(I,state,index,v0)) { n = I->Neighbor[index]; n++; /* skip count */ while(1) { /* look for an attached non-hydrogen as a base */ a1 = I->Neighbor[n]; n+=2; if(a1<0) break; if(ObjectMoleculeGetAtomVertex(I,state,a1,v1)) { subtract3f(v1,v0,n0); normalize3f(n0); /* n0's point away from center atom */ copy3f(n0,occ+3*nOcc); nOcc++; if(nOcc==MaxOcc) /* safety valve */ break; } } switch(ai->geom) { case cAtomInfoPlaner: if(nOcc>1) { cross_product3f(occ,occ+3,normal); if(nOcc>2) { cross_product3f(occ,occ+6,v2); if(dot_product3f(normal,v2)<0) { subtract3f(normal,v2,normal); } else { add3f(normal,v2,normal); } cross_product3f(occ+3,occ+6,v2); if(dot_product3f(normal,v2)<0) { subtract3f(normal,v2,normal); } else { add3f(normal,v2,normal); } } normalize3f(normal); found=true; } break; } } return found; } /*========================================================================*/ int ObjectMoleculeFindOpenValenceVector(ObjectMolecule *I,int state, int index,float *v,float *seek) { CoordSet *cs; int nOcc = 0; float occ[MaxOcc*3]; int last_occ = -1; int n; int a1; float v0[3],v1[3],n0[3],t[3]; int result = false; AtomInfoType *ai,*ai1; float y[3],z[3]; /* default is +X */ v[0]=1.0; v[1]=0.0; v[2]=0.0; if(state<0) state=0; if(I->NCSet==1) state=0; state = state % I->NCSet; cs = I->CSet[state]; if(cs) { if((index>=0)&&(index<=I->NAtom)) { ai=I->AtomInfo+index; if(ObjectMoleculeGetAtomVertex(I,state,index,v0)) { ObjectMoleculeUpdateNeighbors(I); n = I->Neighbor[index]; n++; /* skip count */ while(1) { /* look for an attached non-hydrogen as a base */ a1 = I->Neighbor[n]; n+=2; if(a1<0) break; ai1=I->AtomInfo+a1; if(ObjectMoleculeGetAtomVertex(I,state,a1,v1)) { last_occ = a1; subtract3f(v1,v0,n0); normalize3f(n0); /* n0's point away from center atom */ copy3f(n0,occ+3*nOcc); nOcc++; if(nOcc==MaxOcc) /* safety valve */ break; } } if((!nOcc)||(nOcc>4)||(ai->geom==cAtomInfoNone)) { if(!seek) get_random3f(v); else copy3f(seek,v); result = true; } else { switch(nOcc) { case 1: /* only one current occupied position */ switch(ai->geom) { case cAtomInfoTetrahedral: if(!seek) { get_system1f3f(occ,y,z); scale3f(occ,-0.334F,v); scale3f(z, 0.943F,t); add3f(t,v,v); } else { /* point hydrogen towards sought vector */ copy3f(seek,z); get_system2f3f(occ,z,y); scale3f(occ,-0.334F,v); scale3f(z, 0.943F,t); add3f(t,v,v); } result = true; break; case cAtomInfoPlaner: { if(!seek) { if((last_occ>=0)&&get_planer_normal(I,state,last_occ,n0)) { copy3f(n0,y); get_system2f3f(occ,y,z); } else { get_system1f3f(occ,y,z); } scale3f(occ,-0.500F,v); scale3f(z, 0.866F,t); add3f(t,v,v); } else { copy3f(seek,z); get_system2f3f(occ,z,y); scale3f(occ,-0.500F,v); scale3f(z, 0.866F,t); add3f(t,v,v); } result = true; } break; case cAtomInfoLinear: scale3f(occ,-1.0F,v); result = true; break; default: if(!seek) get_random3f(v); else copy3f(seek,v); result = true; break; } break; case 2: /* only two current occupied positions */ switch(ai->geom) { case cAtomInfoTetrahedral: add3f(occ,occ+3,t); get_system2f3f(t,occ,z); scale3f(t,-1.0F,v); if(seek) { if(dot_product3f(z,seek)<0.0F) { invert3f(z); } } scale3f(z,1.41F,t); add3f(t,v,v); result = true; break; case cAtomInfoPlaner: add3f(occ,occ+3,t); scale3f(t,-1.0F,v); result = true; break; default: if(!seek) get_random3f(v); else copy3f(seek,v); /* hypervalent */ result = true; break; } break; case 3: /* only three current occupied positions */ switch(ai->geom) { case cAtomInfoTetrahedral: add3f(occ,occ+3,t); add3f(occ+6,t,t); scale3f(t,-1.0F,v); result = true; break; default: if(!seek) get_random3f(v); else copy3f(seek,v); /* hypervalent */ result = true; break; } break; case 4: if(!seek) get_random3f(v); else copy3f(seek,v); /* hypervalent */ result = true; break; } } } } } normalize3f(v); return(result); #undef MaxOcc } /*========================================================================*/ void ObjectMoleculeCreateSpheroid(ObjectMolecule *I,int average) { CoordSet *cs; float *spheroid = NULL; int a,b,c,a0; SphereRec *sp; float *spl; float *v,*v0,*s,*f,ang,min_dist,*max_sq; int *i; float *center = NULL; float d0[3],n0[3],d1[3],d2[3]; float p0[3],p1[3],p2[3]; int t0,t1,t2,bt0,bt1,bt2; float dp,l,*fsum = NULL; float *norm = NULL; float spheroid_smooth; float spheroid_fill; float spheroid_ratio=0.1F; /* minimum ratio of width over length */ float spheroid_minimum = 0.02F; /* minimum size - to insure valid normals */ int row,*count=NULL,base; int nRow; int first=0; int last=0; int current; int cscount; int n_state=0; sp=Sphere1; nRow = I->NAtom*sp->nDot; center=Alloc(float,I->NAtom*3); count=Alloc(int,I->NAtom); fsum=Alloc(float,nRow); max_sq = Alloc(float,I->NAtom); spl=spheroid; spheroid_smooth=SettingGet(cSetting_spheroid_smooth); spheroid_fill=SettingGet(cSetting_spheroid_fill); /* first compute average coordinate */ if(average<1) average=I->NCSet; current=0; cscount=0; while(current<I->NCSet) { if(I->CSet[current]) { if(!cscount) first=current; cscount++; last=current+1; } if(cscount==average) { PRINTFB(FB_ObjectMolecule,FB_Details) " ObjectMolecule: computing spheroid from states %d to %d.\n", first+1,last ENDFB; spheroid=Alloc(float,nRow); v=center; i = count; for(a=0;a<I->NAtom;a++) { *(v++)=0.0; *(v++)=0.0; *(v++)=0.0; *(i++)=0; } for(b=first;b<last;b++) { cs=I->CSet[b]; if(cs) { v = cs->Coord; for(a=0;a<cs->NIndex;a++) { a0=cs->IdxToAtm[a]; v0 = center+3*a0; add3f(v,v0,v0); (*(count+a0))++; v+=3; } } } i=count; v=center; for(a=0;a<I->NAtom;a++) if(*i) { (*(v++))/=(*i); (*(v++))/=(*i); (*(v++))/=(*i++); } else { v+=3; i++; } /* now go through and compute radial distances */ f = fsum; s = spheroid; for(a=0;a<nRow;a++) { *(f++)=0.0; *(s++)=0.0; } v = max_sq; for(a=0;a<I->NAtom;a++) *(v++)=0.0; for(b=first;b<last;b++) { cs=I->CSet[b]; if(cs) { v = cs->Coord; for(a=0;a<cs->NIndex;a++) { a0=cs->IdxToAtm[a]; base = (a0*sp->nDot); v0 = center+(3*a0); subtract3f(v,v0,d0); /* subtract from average */ l = lengthsq3f(d0); if(l>max_sq[a0]) max_sq[a0]=l; if(l>0.0) { float isq = (float)(1.0/sqrt1d(l)); scale3f(d0,isq,n0); for(c=0;c<sp->nDot;c++) { /* average over spokes */ dp=dot_product3f(sp->dot[c],n0); row = base + c; if(dp>=0.0) { ang = (float)((acos(dp)/spheroid_smooth)*(cPI/2.0)); if(ang>spheroid_fill) ang=spheroid_fill; /* take envelop to zero over that angle */ if(ang<=(cPI/2.0)) { dp = (float)cos(ang); fsum[row] += dp*dp; spheroid[row] += l*dp*dp*dp; } } } } v+=3; } } } f=fsum; s=spheroid; for(a=0;a<I->NAtom;a++) { min_dist = (float)(spheroid_ratio*sqrt(max_sq[a])); if(min_dist<spheroid_minimum) min_dist=spheroid_minimum; for(b=0;b<sp->nDot;b++) { if(*f>R_SMALL4) { (*s)=(float)(sqrt1d((*s)/(*(f++)))); /* we put the "rm" in "rms" */ } else { f++; } if(*s<min_dist) *s=min_dist; s++; } } /* set frame 0 coordinates to the average */ cs=I->CSet[first]; if(cs) { v = cs->Coord; for(a=0;a<cs->NIndex;a++) { a0=cs->IdxToAtm[a]; v0 = center+3*a0; copy3f(v0,v); v+=3; } } /* now compute surface normals */ norm = Alloc(float,nRow*3); for(a=0;a<nRow;a++) { zero3f(norm+a*3); } for(a=0;a<I->NAtom;a++) { base = a*sp->nDot; for(b=0;b<sp->NTri;b++) { t0 = sp->Tri[b*3 ]; t1 = sp->Tri[b*3+1]; t2 = sp->Tri[b*3+2]; bt0 = base + t0; bt1 = base + t1; bt2 = base + t2; copy3f(sp->dot[t0],p0); copy3f(sp->dot[t1],p1); copy3f(sp->dot[t2],p2); /* scale3f(sp->dot[t0].v,spheroid[bt0],p0); scale3f(sp->dot[t1].v,spheroid[bt1],p1); scale3f(sp->dot[t2].v,spheroid[bt2],p2);*/ subtract3f(p1,p0,d1); subtract3f(p2,p0,d2); cross_product3f(d1,d2,n0); normalize3f(n0); v = norm+bt0*3; add3f(n0,v,v); v = norm+bt1*3; add3f(n0,v,v); v = norm+bt2*3; add3f(n0,v,v); } } f=norm; for(a=0;a<I->NAtom;a++) { base = a*sp->nDot; for(b=0;b<sp->nDot;b++) { normalize3f(f); f+=3; } } if(I->CSet[first]) { I->CSet[first]->Spheroid=spheroid; I->CSet[first]->SpheroidNormal=norm; I->CSet[first]->NSpheroid=nRow; } else { FreeP(spheroid); FreeP(norm); } for(b=first+1;b<last;b++) { cs=I->CSet[b]; if(cs) { if(cs->fFree) cs->fFree(cs); } I->CSet[b]=NULL; } if(n_state!=first) { I->CSet[n_state]=I->CSet[first]; I->CSet[first]=NULL; } n_state++; cscount=0; } current++; } I->NCSet=n_state; FreeP(center); FreeP(count); FreeP(fsum); FreeP(max_sq); ObjectMoleculeInvalidate(I,cRepSphere,cRepInvProp); } /*========================================================================*/ void ObjectMoleculeReplaceAtom(ObjectMolecule *I,int index,AtomInfoType *ai) { if((index>=0)&&(index<=I->NAtom)) { memcpy(I->AtomInfo+index,ai,sizeof(AtomInfoType)); ObjectMoleculeInvalidate(I,cRepAll,cRepInvAtoms); /* could we put in a refinement step here? */ } } /*========================================================================*/ void ObjectMoleculePrepareAtom(ObjectMolecule *I,int index,AtomInfoType *ai) { /* match existing properties of the old atom */ int a; AtomInfoType *ai0; if((index>=0)&&(index<=I->NAtom)) { ai0=I->AtomInfo + index; ai->resv=ai0->resv; ai->hetatm=ai0->hetatm; ai->flags=ai0->flags; ai->geom=ai0->geom; /* ?*/ ai->q=ai0->q; ai->b=ai0->b; strcpy(ai->chain,ai0->chain); strcpy(ai->alt,ai0->alt); strcpy(ai->resi,ai0->resi); strcpy(ai->segi,ai0->segi); strcpy(ai->resn,ai0->resn); if((ai->elem[0]==ai0->elem[0])&&(ai->elem[1]==ai0->elem[1])) ai->color=ai0->color; else if((ai->elem[0]=='C')&&(ai->elem[1]==0)) /* carbons are always colored according to the object color */ ai->color=I->Obj.Color; else ai->color=AtomInfoGetColor(ai); for(a=0;a<cRepCnt;a++) ai->visRep[a]=ai0->visRep[a]; ai->id=-1; AtomInfoUniquefyNames(I->AtomInfo,I->NAtom,ai,1); AtomInfoAssignParameters(ai); } } /*========================================================================*/ void ObjectMoleculePreposReplAtom(ObjectMolecule *I,int index, AtomInfoType *ai) { int n; int a1; AtomInfoType *ai1; float v0[3],v1[3],v[3]; float d0[3],d,n0[3]; int cnt; float t[3],sum[3]; int a; int ncycle; ObjectMoleculeUpdateNeighbors(I); for(a=0;a<I->NCSet;a++) { if(I->CSet[a]) { if(ObjectMoleculeGetAtomVertex(I,a,index,v0)) { copy3f(v0,v); /* default is direct superposition */ ncycle=-1; while(ncycle) { cnt = 0; n = I->Neighbor[index]; n++; /* skip count */ zero3f(sum); while(1) { /* look for an attached non-hydrogen as a base */ a1 = I->Neighbor[n]; n+=2; if(a1<0) break; ai1=I->AtomInfo+a1; if(ai1->protons!=1) if(ObjectMoleculeGetAtomVertex(I,a,a1,v1)) { d = AtomInfoGetBondLength(ai,ai1); subtract3f(v0,v1,n0); normalize3f(n0); scale3f(n0,d,d0); add3f(d0,v1,t); add3f(t,sum,sum); cnt++; } } if(cnt) { scale3f(sum,1.0F/cnt,sum); copy3f(sum,v0); if((cnt>1)&&(ncycle<0)) ncycle=5; } ncycle=abs(ncycle)-1; } if(cnt) copy3f(sum,v); ObjectMoleculeSetAtomVertex(I,a,index,v); } } } } /*========================================================================*/ void ObjectMoleculeSaveUndo(ObjectMolecule *I,int state,int log) { CoordSet *cs; FreeP(I->UndoCoord[I->UndoIter]); I->UndoState[I->UndoIter]=-1; if(state<0) state=0; if(I->NCSet==1) state=0; state = state % I->NCSet; cs = I->CSet[state]; if(cs) { I->UndoCoord[I->UndoIter] = Alloc(float,cs->NIndex*3); memcpy(I->UndoCoord[I->UndoIter],cs->Coord,sizeof(float)*cs->NIndex*3); I->UndoState[I->UndoIter]=state; I->UndoNIndex[I->UndoIter] = cs->NIndex; } I->UndoIter=cUndoMask&(I->UndoIter+1); ExecutiveSetLastObjectEdited((CObject*)I); if(log) { OrthoLineType line; if(SettingGet(cSetting_logging)) { sprintf(line,"cmd.push_undo(\"%s\",%d)\n",I->Obj.Name,state+1); PLog(line,cPLog_no_flush); } } } /*========================================================================*/ void ObjectMoleculeUndo(ObjectMolecule *I,int dir) { CoordSet *cs; int state; FreeP(I->UndoCoord[I->UndoIter]); I->UndoState[I->UndoIter]=-1; state=SceneGetState(); if(state<0) state=0; if(I->NCSet==1) state=0; state = state % I->NCSet; cs = I->CSet[state]; if(cs) { I->UndoCoord[I->UndoIter] = Alloc(float,cs->NIndex*3); memcpy(I->UndoCoord[I->UndoIter],cs->Coord,sizeof(float)*cs->NIndex*3); I->UndoState[I->UndoIter]=state; I->UndoNIndex[I->UndoIter] = cs->NIndex; } I->UndoIter=cUndoMask&(I->UndoIter+dir); if(!I->UndoCoord[I->UndoIter]) I->UndoIter=cUndoMask&(I->UndoIter-dir); if(I->UndoState[I->UndoIter]>=0) { state=I->UndoState[I->UndoIter]; if(state<0) state=0; if(I->NCSet==1) state=0; state = state % I->NCSet; cs = I->CSet[state]; if(cs) { if(cs->NIndex==I->UndoNIndex[I->UndoIter]) { memcpy(cs->Coord,I->UndoCoord[I->UndoIter],sizeof(float)*cs->NIndex*3); I->UndoState[I->UndoIter]=-1; FreeP(I->UndoCoord[I->UndoIter]); if(cs->fInvalidateRep) cs->fInvalidateRep(cs,cRepAll,cRepInvCoord); SceneChanged(); } } } } /*========================================================================*/ int ObjectMoleculeAddBond(ObjectMolecule *I,int sele0,int sele1,int order) { int a1,a2; AtomInfoType *ai1,*ai2; int s1,s2; int c = 0; BondType *bnd; ai1=I->AtomInfo; for(a1=0;a1<I->NAtom;a1++) { s1=ai1->selEntry; if(SelectorIsMember(s1,sele0)) { ai2=I->AtomInfo; for(a2=0;a2<I->NAtom;a2++) { s2=ai2->selEntry; if(SelectorIsMember(s2,sele1)) { { VLACheck(I->Bond,BondType,I->NBond); bnd = I->Bond+(I->NBond); bnd->index[0]=a1; bnd->index[1]=a2; bnd->order=order; bnd->stereo = 0; bnd->id=-1; I->NBond++; c++; I->AtomInfo[a1].chemFlag=false; I->AtomInfo[a2].chemFlag=false; } } ai2++; } } ai1++; } if(c) { ObjectMoleculeInvalidate(I,cRepLine,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepCyl,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepNonbonded,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepNonbondedSphere,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepRibbon,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepCartoon,cRepInvBonds); ObjectMoleculeUpdateIDNumbers(I); } return(c); } /*========================================================================*/ int ObjectMoleculeAdjustBonds(ObjectMolecule *I,int sele0,int sele1,int mode,int order) { int a0,a1; int offset=0; BondType *b0; int both; int s; int a; offset=0; b0=I->Bond; for(a=0;a<I->NBond;a++) { a0=b0->index[0]; a1=b0->index[1]; both=0; s=I->AtomInfo[a0].selEntry; if(SelectorIsMember(s,sele0)) both++; s=I->AtomInfo[a1].selEntry; if(SelectorIsMember(s,sele1)) both++; if(both<2) { /* reverse combo */ both=0; s=I->AtomInfo[a1].selEntry; if(SelectorIsMember(s,sele0)) both++; s=I->AtomInfo[a0].selEntry; if(SelectorIsMember(s,sele1)) both++; } if(both==2) { switch(mode) { case 0: /* cycle */ b0->order++; if(b0->order>3) b0->order=1; I->AtomInfo[a0].chemFlag=false; I->AtomInfo[a1].chemFlag=false; break; case 1: /* set */ b0->order=order; I->AtomInfo[a0].chemFlag=false; I->AtomInfo[a1].chemFlag=false; break; } ObjectMoleculeInvalidate(I,cRepLine,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepCyl,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepNonbonded,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepNonbondedSphere,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepRibbon,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepCartoon,cRepInvBonds); } b0++; } return(-offset); } /*========================================================================*/ int ObjectMoleculeRemoveBonds(ObjectMolecule *I,int sele0,int sele1) { int a0,a1; int offset=0; BondType *b0,*b1; int both; int s; int a; offset=0; b0=I->Bond; b1=I->Bond; for(a=0;a<I->NBond;a++) { a0=b0->index[0]; a1=b0->index[1]; both=0; s=I->AtomInfo[a0].selEntry; if(SelectorIsMember(s,sele0)) both++; s=I->AtomInfo[a1].selEntry; if(SelectorIsMember(s,sele1)) both++; if(both<2) { /* reverse combo */ both=0; s=I->AtomInfo[a1].selEntry; if(SelectorIsMember(s,sele0)) both++; s=I->AtomInfo[a0].selEntry; if(SelectorIsMember(s,sele1)) both++; } if(both==2) { offset--; b0++; I->AtomInfo[a0].chemFlag=false; I->AtomInfo[a1].chemFlag=false; } else if(offset) { *(b1++)=*(b0++); /* copy bond info */ } else { *(b1++)=*(b0++); /* copy bond info */ } } if(offset) { I->NBond += offset; VLASize(I->Bond,BondType,I->NBond); ObjectMoleculeInvalidate(I,cRepLine,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepCyl,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepNonbonded,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepNonbondedSphere,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepRibbon,cRepInvBonds); ObjectMoleculeInvalidate(I,cRepCartoon,cRepInvBonds); } return(-offset); } /*========================================================================*/ void ObjectMoleculePurge(ObjectMolecule *I) { int a,a0,a1; int *oldToNew = NULL; int offset=0; BondType *b0,*b1; AtomInfoType *ai0,*ai1; PRINTFD(FB_ObjectMolecule) " ObjMolPurge-Debug: step 1, delete object selection\n" ENDFD; SelectorDelete(I->Obj.Name); /* remove the object selection and free up any selection entries*/ /* note that we don't delete atom selection members -- those may be needed in the new object */ PRINTFD(FB_ObjectMolecule) " ObjMolPurge-Debug: step 2, purge coordinate sets\n" ENDFD; for(a=0;a<I->NCSet;a++) if(I->CSet[a]) CoordSetPurge(I->CSet[a]); if(I->CSTmpl) { CoordSetPurge(I->CSTmpl); } PRINTFD(FB_ObjectMolecule) " ObjMolPurge-Debug: step 3, old-to-new mapping\n" ENDFD; oldToNew = Alloc(int,I->NAtom); ai0=I->AtomInfo; ai1=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(ai0->deleteFlag) { offset--; ai0++; oldToNew[a]=-1; } else if(offset) { *(ai1++)=*(ai0++); oldToNew[a]=a+offset; } else { oldToNew[a]=a; ai0++; ai1++; } } if(offset) { I->NAtom += offset; VLASize(I->AtomInfo,AtomInfoType,I->NAtom); for(a=0;a<I->NCSet;a++) if(I->CSet[a]) CoordSetAdjustAtmIdx(I->CSet[a],oldToNew,I->NAtom); } PRINTFD(FB_ObjectMolecule) " ObjMolPurge-Debug: step 4, bonds\n" ENDFD; offset=0; b0=I->Bond; b1=I->Bond; for(a=0;a<I->NBond;a++) { a0=b0->index[0]; a1=b0->index[1]; if((oldToNew[a0]<0)||(oldToNew[a1]<0)) { offset--; b0++; } else if(offset) { *b1=*b0; b1->index[0]=oldToNew[a0]; /* copy bond info */ b1->index[1]=oldToNew[a1]; b0++; b1++; } else { *b1=*b0; b1->index[0]=oldToNew[a0]; /* copy bond info */ b1->index[1]=oldToNew[a1]; b0++; b1++; /* TODO check reasoning agaist above */ } } if(offset) { I->NBond += offset; VLASize(I->Bond,BondType,I->NBond); } FreeP(oldToNew); PRINTFD(FB_ObjectMolecule) " ObjMolPurge-Debug: step 5, invalidate...\n" ENDFD; ObjectMoleculeInvalidate(I,cRepAll,cRepInvAtoms); PRINTFD(FB_ObjectMolecule) " ObjMolPurge-Debug: leaving...\n" ENDFD; } /*========================================================================*/ int ObjectMoleculeGetAtomGeometry(ObjectMolecule *I,int state,int at) { /* this determines hybridization from coordinates in those few cases * where it is unambiguous */ int result = -1; int n,nn; float v0[3],v1[3],v2[3],v3[3]; float d1[3],d2[3],d3[3]; float cp1[3],cp2[3],cp3[3]; float avg; float dp; n = I->Neighbor[at]; nn = I->Neighbor[n++]; /* get count */ if(nn==4) result = cAtomInfoTetrahedral; else if(nn==3) { /* check cross products */ ObjectMoleculeGetAtomVertex(I,state,at,v0); ObjectMoleculeGetAtomVertex(I,state,I->Neighbor[n],v1); ObjectMoleculeGetAtomVertex(I,state,I->Neighbor[n+2],v2); ObjectMoleculeGetAtomVertex(I,state,I->Neighbor[n+4],v3); subtract3f(v1,v0,d1); subtract3f(v2,v0,d2); subtract3f(v3,v0,d3); cross_product3f(d1,d2,cp1); cross_product3f(d2,d3,cp2); cross_product3f(d3,d1,cp3); normalize3f(cp1); normalize3f(cp2); normalize3f(cp3); avg=(dot_product3f(cp1,cp2)+ dot_product3f(cp2,cp3)+ dot_product3f(cp3,cp1))/3.0F; if(avg>0.75) result=cAtomInfoPlaner; else result=cAtomInfoTetrahedral; } else if(nn==2) { ObjectMoleculeGetAtomVertex(I,state,at,v0); ObjectMoleculeGetAtomVertex(I,state,I->Neighbor[n],v1); ObjectMoleculeGetAtomVertex(I,state,I->Neighbor[n+2],v2); subtract3f(v1,v0,d1); subtract3f(v2,v0,d2); normalize3f(d1); normalize3f(d2); dp = dot_product3f(d1,d2); if(dp<-0.75) result=cAtomInfoLinear; } return(result); } /*========================================================================*/ void ObjectMoleculeInferChemForProtein(ObjectMolecule *I,int state) { /* Infers chemical relations for a molecules under protein assumptions. * * NOTE: this routine needs an all-atom model (with hydrogens!) * and it will make mistakes on non-protein atoms (if they haven't * already been assigned) */ int a,n,a0,a1,nn; int changedFlag = true; AtomInfoType *ai,*ai0,*ai1=NULL; ObjectMoleculeUpdateNeighbors(I); /* first, try to find all amids and acids */ while(changedFlag) { changedFlag=false; for(a=0;a<I->NAtom;a++) { ai=I->AtomInfo+a; if(ai->chemFlag) { if(ai->geom==cAtomInfoPlaner) if(ai->protons == cAN_C) { n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(nn>1) { a1 = -1; while(1) { a0 = I->Neighbor[n]; n+=2; if(a0<0) break; ai0 = I->AtomInfo+a0; if((ai0->protons==cAN_O)&&(!ai0->chemFlag)) { a1=a0; ai1=ai0; /* found candidate carbonyl */ break; } } if(a1>0) { n = I->Neighbor[a]+1; while(1) { a0 = I->Neighbor[n]; if(a0<0) break; n+=2; if(a0!=a1) { ai0 = I->AtomInfo+a0; if(ai0->protons==cAN_O) { if(!ai0->chemFlag) { ai0->chemFlag=true; /* acid */ ai0->geom=cAtomInfoPlaner; ai0->valence=1; ai1->chemFlag=true; ai1->geom=cAtomInfoPlaner; ai1->valence=1; changedFlag=true; break; } } else if(ai0->protons==cAN_N) { if(!ai0->chemFlag) { ai0->chemFlag=true; /* amide N */ ai0->geom=cAtomInfoPlaner; ai0->valence=3; ai1->chemFlag=true; /* amide =O */ ai1->geom=cAtomInfoPlaner; ai1->valence=1; changedFlag=true; break; } else if(ai0->geom==cAtomInfoPlaner) { ai1->chemFlag=true; /* amide =O */ ai1->geom=cAtomInfoPlaner; ai1->valence=1; changedFlag=true; break; } } } } } } } } } } /* then handle aldehydes and amines (partial amides - both missing a valence) */ changedFlag=true; while(changedFlag) { changedFlag=false; for(a=0;a<I->NAtom;a++) { ai=I->AtomInfo+a; if(!ai->chemFlag) { if(ai->protons==cAN_C) { n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(nn>1) { a1 = -1; while(1) { a0 = I->Neighbor[n]; n+=2; if(a0<0) break; ai0 = I->AtomInfo+a0; if((ai0->protons==cAN_O)&&(!ai0->chemFlag)) { /* =O */ ai->chemFlag=true; ai->geom=cAtomInfoPlaner; ai->valence=1; ai0->chemFlag=true; ai0->geom=cAtomInfoPlaner; ai0->valence=3; changedFlag=true; break; } } } } else if(ai->protons==cAN_N) { if((!ai->chemFlag)||ai->geom!=cAtomInfoLinear) { if(ai->formalCharge==0.0) { ai->chemFlag=true; ai->geom=cAtomInfoPlaner; ai->valence=3; } } } } } } } /*========================================================================*/ void ObjectMoleculeInferChemFromNeighGeom(ObjectMolecule *I,int state) { /* infers chemical relations from neighbors and geometry * NOTE: very limited in scope */ int a,n,a0,nn; int changedFlag=true; int geom; int carbonVal[10]; AtomInfoType *ai,*ai2; carbonVal[cAtomInfoTetrahedral] = 4; carbonVal[cAtomInfoPlaner] = 3; carbonVal[cAtomInfoLinear] = 2; ObjectMoleculeUpdateNeighbors(I); while(changedFlag) { changedFlag=false; for(a=0;a<I->NAtom;a++) { ai=I->AtomInfo+a; if(!ai->chemFlag) { geom=ObjectMoleculeGetAtomGeometry(I,state,a); switch(ai->protons) { case cAN_K: ai->chemFlag=1; ai->geom=cAtomInfoNone; ai->valence=0; break; case cAN_H: case cAN_F: case cAN_I: case cAN_Br: ai->chemFlag=1; ai->geom=cAtomInfoSingle; ai->valence=1; break; case cAN_O: n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(nn!=1) { /* water, hydroxy, ether */ ai->chemFlag=1; ai->geom=cAtomInfoTetrahedral; ai->valence=2; } else { /* hydroxy or carbonyl? check carbon geometry */ a0 = I->Neighbor[n+2]; ai2=I->AtomInfo+a0; if(ai2->chemFlag) { if((ai2->geom==cAtomInfoTetrahedral)|| (ai2->geom==cAtomInfoLinear)) { ai->chemFlag=1; /* hydroxy */ ai->geom=cAtomInfoTetrahedral; ai->valence=2; } } } break; case cAN_C: if(geom>=0) { ai->geom = geom; ai->valence = carbonVal[geom]; ai->chemFlag=true; } else { n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(nn==1) { /* only one neighbor */ ai2=I->AtomInfo+I->Neighbor[n]; if(ai2->chemFlag&&(ai2->geom==cAtomInfoTetrahedral)) { ai->chemFlag=true; /* singleton carbon bonded to tetC must be tetC */ ai->geom=cAtomInfoTetrahedral; ai->valence=4; } } } break; case cAN_N: if(geom==cAtomInfoPlaner) { ai->chemFlag=true; ai->geom=cAtomInfoPlaner; ai->valence=3; } else if(geom==cAtomInfoTetrahedral) { ai->chemFlag=true; ai->geom=cAtomInfoTetrahedral; ai->valence=4; } break; case cAN_S: n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(nn==4) { /* sulfone */ ai->chemFlag=true; ai->geom=cAtomInfoTetrahedral; ai->valence=4; } else if(nn==3) { /* suloxide */ ai->chemFlag=true; ai->geom=cAtomInfoTetrahedral; ai->valence=3; } else if(nn==2) { /* thioether */ ai->chemFlag=true; ai->geom=cAtomInfoTetrahedral; ai->valence=2; } break; case cAN_Cl: ai->chemFlag=1; if(ai->formalCharge==0.0) { ai->geom=cAtomInfoSingle; ai->valence=1; } else { ai->geom=cAtomInfoNone; ai->valence=0; } break; } if(ai->chemFlag) changedFlag=true; } } } } /*========================================================================*/ void ObjectMoleculeInferHBondFromChem(ObjectMolecule *I) { int a; AtomInfoType *ai; int a1; int n,nn; int has_hydro,may_have_lone_pair; /* initialize accumulators on uncategorized atoms */ ObjectMoleculeUpdateNeighbors(I); ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { n = I->Neighbor[a]; nn = I->Neighbor[n++]; ai->hb_donor = false; ai->hb_acceptor = false; has_hydro = (nn < ai->valence); /* implicit hydrogens? */ if(!has_hydro) { /* explicit hydrogens? */ has_hydro = false; switch(ai->protons) { case cAN_N: case cAN_O: while((a1 = I->Neighbor[n])>=0) { n+=2; if(I->AtomInfo[a1].protons==1) { has_hydro=true; break; } } break; } } switch(ai->protons) { /* cat-ions, lewis acids etc. */ case cAN_Fe: case cAN_Ca: case cAN_Cu: case cAN_K: case cAN_Na: case cAN_Mg: case cAN_Zn: case cAN_Hg: ai->hb_donor=true; break; case cAN_N: if(has_hydro) ai->hb_donor=true; else { may_have_lone_pair = false; if((!has_hydro)&&(ai->protons==cAN_N)) { n = I->Neighbor[a] + 1; while(I->Neighbor[n]>=0) { if(I->Neighbor[n+1]>1) { /* any double/triple/delocalized bonds? */ may_have_lone_pair = true; } n+=2; } } if((ai->formalCharge<=0)&&may_have_lone_pair) { ai->hb_acceptor = true; } } break; case cAN_O: if(has_hydro) ai->hb_donor = true; if(ai->formalCharge<=0) ai->hb_acceptor = true; break; } ai++; } } /*========================================================================*/ void ObjectMoleculeInferChemFromBonds(ObjectMolecule *I,int state) { int a,b; BondType *b0; AtomInfoType *ai,*ai0,*ai1=NULL; int a0,a1; int expect,order; int n,nn; int changedFlag; /* initialize accumulators on uncategorized atoms */ ObjectMoleculeUpdateNeighbors(I); ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(!ai->chemFlag) { ai->geom=0; ai->valence=0; } ai++; } /* find maximum bond order for each atom */ b0=I->Bond; for(b=0;b<I->NBond;b++) { a0 = b0->index[0]; a1 = b0->index[1]; ai0=I->AtomInfo + a0; ai1=I->AtomInfo + a1; order = b0->order; b0++; if(!ai0->chemFlag) { if(order>ai0->geom) ai0->geom=order; ai0->valence+=order; } if(!ai1->chemFlag) { if(order>ai1->geom) ai1->geom=order; ai1->valence+=order; } if(order==3) { /* override existing chemistry * this is a temp fix to a pressing problem... we need to rethink the chemisty assignment ordering (should bond information come first? */ ai0->geom = cAtomInfoLinear; ai1->geom = cAtomInfoLinear; switch(ai0->protons) { case cAN_C: ai0->valence=2; break; default: ai0->valence=1; } switch(ai1->protons) { case cAN_C: ai1->valence=2; break; default: ai1->valence=1; } ai0->chemFlag=true; ai1->chemFlag=true; } else if(order==4) { ai0->geom = cAtomInfoPlaner; ai1->geom = cAtomInfoPlaner; switch(ai0->protons) { case cAN_O: ai0->valence=1; break; case cAN_N: ai0->valence=2; break; case cAN_C: ai0->valence=3; break; default: ai0->valence=4; } switch(ai1->protons) { case cAN_O: ai1->valence=1; break; case cAN_N: ai1->valence=2; break; case cAN_C: ai1->valence=3; break; default: ai1->valence=1; } ai0->chemFlag=true; ai1->chemFlag=true; } } /* now set up valences and geometries */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(!ai->chemFlag) { expect = AtomInfoGetExpectedValence(ai); n = I->Neighbor[a]; nn = I->Neighbor[n++]; if(ai->geom==3) { ai->geom = cAtomInfoLinear; switch(ai->protons) { case cAN_C: ai->valence=2; break; default: ai->valence=1; } ai->chemFlag=true; } else { if(expect<0) expect = -expect; /* for now, just ignore this issue */ /* printf("%d %d %d %d\n",ai->geom,ai->valence,nn,expect);*/ if(ai->valence==expect) { /* sum of bond orders equals valence */ ai->chemFlag=true; ai->valence=nn; switch(ai->geom) { /* max bond order observed */ case 0: ai->geom = cAtomInfoNone; break; case 2: ai->geom = cAtomInfoPlaner; break; case 3: ai->geom = cAtomInfoLinear; break; default: if(expect==1) ai->geom = cAtomInfoSingle; else ai->geom = cAtomInfoTetrahedral; break; } } else if(ai->valence<expect) { /* missing a bond */ ai->chemFlag=true; ai->valence=nn+(expect-ai->valence); switch(ai->geom) { case 2: ai->geom = cAtomInfoPlaner; break; case 3: ai->geom = cAtomInfoLinear; break; default: if(expect==1) ai->geom = cAtomInfoSingle; else ai->geom = cAtomInfoTetrahedral; break; } } else if(ai->valence>expect) { ai->chemFlag=true; ai->valence=nn; switch(ai->geom) { case 2: ai->geom = cAtomInfoPlaner; break; case 3: ai->geom = cAtomInfoLinear; break; default: if(expect==1) ai->geom = cAtomInfoSingle; else ai->geom = cAtomInfoTetrahedral; break; } if(nn>3) ai->geom = cAtomInfoTetrahedral; } } } ai++; } /* now go through and make sure conjugated amines are planer */ changedFlag=true; while(changedFlag) { changedFlag=false; ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(ai->chemFlag) { if(ai->protons==cAN_N) if(ai->formalCharge==0) if(ai->geom==cAtomInfoTetrahedral) { /* search for uncharged tetrahedral nitrogen */ n = I->Neighbor[a]+1; while(1) { a0 = I->Neighbor[n]; n+=2; if(a0<0) break; ai0 = I->AtomInfo+a0; if((ai0->chemFlag)&&(ai0->geom==cAtomInfoPlaner)&& ((ai0->protons==cAN_C)||(ai0->protons==cAN_N))) { ai->geom=cAtomInfoPlaner; /* found probable delocalization */ ai->valence=3; /* just in case...*/ changedFlag=true; break; } } } } ai++; } } /* now go through and make sure conjugated anions are planer */ changedFlag=true; while(changedFlag) { changedFlag=false; ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(ai->chemFlag) { if(ai->protons==cAN_O) if(ai->formalCharge==-1) if((ai->geom==cAtomInfoTetrahedral)|| (ai->geom==cAtomInfoSingle)) { /* search for anionic tetrahedral oxygen */ n = I->Neighbor[a]+1; while(1) { a0 = I->Neighbor[n]; n+=2; if(a0<0) break; ai0 = I->AtomInfo+a0; if((ai0->chemFlag)&&(ai0->geom==cAtomInfoPlaner)&& ((ai0->protons==cAN_C)||(ai0->protons==cAN_N))) { ai->geom=cAtomInfoPlaner; /* found probable delocalization */ changedFlag=true; break; } } } } ai++; } } } /*========================================================================*/ int ObjectMoleculeTransformSelection(ObjectMolecule *I,int state, int sele,float *TTT,int log,char *sname) { /* if sele == -1, then the whole object state is transformed */ int a,s; int flag=false; CoordSet *cs; AtomInfoType *ai; int logging; int all_states=false,inp_state; int ok=true; inp_state=state; if(state==-1) state=ObjectGetCurrentState(&I->Obj,false); if(state<0) { all_states=true; state=-1; } PRINTFD(FB_ObjectMolecule) "ObjMolTransSele-Debug: state %d\n",state ENDFD; while(1) { if(all_states) { state++; if(state>=I->NCSet) break; } if(state<I->NCSet) { cs = I->CSet[state]; if(cs) { if(sele>=0) { ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(!(ai->protekted==1)) if(SelectorIsMember(s,sele)) { CoordSetTransformAtom(cs,a,TTT); flag=true; } ai++; } } else { ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { if(!(ai->protekted==1)) CoordSetTransformAtom(cs,a,TTT); ai++; } flag=true; } if(flag) cs->fInvalidateRep(cs,cRepAll,cRepInvCoord); } } if(!all_states) break; } if(log) { OrthoLineType line; WordType sele_str = ",'"; logging = (int)SettingGet(cSetting_logging); if(sele>=0) { strcat(sele_str,sname); strcat(sele_str,"'"); } else sele_str[0]=0; switch(logging) { case cPLog_pml: sprintf(line, "_ cmd.transform_object('%s',[\\\n_ %15.9f,%15.9f,%15.9f,%15.9f,\\\n_ %15.9f,%15.9f,%15.9f,%15.9f,\\\n_ %15.9f,%15.9f,%15.9f,%15.9f,\\\n_ %15.9f,%15.9f,%15.9f,%15.9f\\\n_ ],%d,%d%s)\n", I->Obj.Name, TTT[ 0],TTT[ 1],TTT[ 2],TTT[ 3], TTT[ 4],TTT[ 5],TTT[ 6],TTT[ 7], TTT[ 8],TTT[ 9],TTT[10],TTT[11], TTT[12],TTT[13],TTT[14],TTT[15],inp_state+1,log,sele_str); PLog(line,cPLog_no_flush); break; case cPLog_pym: sprintf(line, "cmd.transform_object('%s',[\n%15.9f,%15.9f,%15.9f,%15.9f,\n%15.9f,%15.9f,%15.9f,%15.9f,\n%15.9f,%15.9f,%15.9f,%15.9f,\n%15.9f,%15.9f,%15.9f,%15.9f\n],%d,%d%s)\n", I->Obj.Name, TTT[ 0],TTT[ 1],TTT[ 2],TTT[ 3], TTT[ 4],TTT[ 5],TTT[ 6],TTT[ 7], TTT[ 8],TTT[ 9],TTT[10],TTT[11], TTT[12],TTT[13],TTT[14],TTT[15],inp_state+1,log,sele_str); PLog(line,cPLog_no_flush); break; default: break; } } return(ok); } /*========================================================================*/ int ObjectMoleculeGetAtomIndex(ObjectMolecule *I,int sele) { int a,s; if(sele<0) return(-1); for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) return(a); } return(-1); } /*========================================================================*/ void ObjectMoleculeUpdateNonbonded(ObjectMolecule *I) { int a; BondType *b; AtomInfoType *ai; /* if(!I->DiscreteFlag) { */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) (ai++)->bonded = false; b=I->Bond; ai=I->AtomInfo; for(a=0;a<I->NBond;a++) { ai[b->index[0]].bonded=true; ai[b->index[1]].bonded=true; b++; } /* } */ } /*========================================================================*/ int ObjectMoleculeGetTotalAtomValence(ObjectMolecule *I,int atom) { int result = 0; int n0; ObjectMoleculeUpdateNeighbors(I); if(atom<I->NAtom) { n0 = I->Neighbor[atom]+1; while(I->Neighbor[n0]>=0) { result+=I->Neighbor[n0+1]; } n0+=2; } else { result = -1; /* error */ } return result; } /*========================================================================*/ void ObjectMoleculeUpdateNeighbors(ObjectMolecule *I) { /* neighbor storage structure: VERY COMPLICATED... 0 list offset for atom 0 = n 1 list offset for atom 1 = n + m + 1 ... n-1 list offset for atom n-1 n count for atom 0 n+1 neighbor of atom 0 n+2 bond index n+3 neighbor of atom 0 n+4 bond index ... n+m -1 terminator for atom 0 n+m+1 count for atom 1 n+m+2 neighbor of atom 1 n+m+3 bond index n+m+4 neighbor of atom 1 n+m+5 bond index etc. NOTE: all atoms have an offset and a terminator whether or not they have any bonds */ int size; int a,b,c,d,l0,l1,*l; BondType *bnd; if(!I->Neighbor) { size = (I->NAtom*3)+(I->NBond*4); if(I->Neighbor) { VLACheck(I->Neighbor,int,size); } else { I->Neighbor=VLAlloc(int,size); } /* initialize */ l = I->Neighbor; for(a=0;a<I->NAtom;a++) (*l++)=0; /* count neighbors for each atom */ bnd = I->Bond; for(b=0;b<I->NBond;b++) { I->Neighbor[bnd->index[0]]++; I->Neighbor[bnd->index[1]]++; bnd++; } /* set up offsets and list terminators */ c = I-> NAtom; for(a=0;a<I->NAtom;a++) { d = I->Neighbor[a]; /* get number of neighbors */ I->Neighbor[c]=d; /* store neighbor count */ I->Neighbor[a]=c+d+d+1; /* set initial position to end of list, we'll fill backwards */ I->Neighbor[I->Neighbor[a]]=-1; /* store terminator */ c += d + d + 2; } /* now load neighbors in a sequential list for each atom (reverse order) */ bnd = I->Bond; for(b=0;b<I->NBond;b++) { l0 = bnd->index[0]; l1 = bnd->index[1]; bnd++; I->Neighbor[l0]--; I->Neighbor[I->Neighbor[l0]]=b; /* store bond indices (for I->Bond) */ I->Neighbor[l0]--; I->Neighbor[I->Neighbor[l0]]=l1; /* store neighbor references (I->AtomInfo, etc.)*/ I->Neighbor[l1]--; I->Neighbor[I->Neighbor[l1]]=b; /* store bond indices (for I->Bond) */ I->Neighbor[l1]--; I->Neighbor[I->Neighbor[l1]]=l0; /* store neighbor references (I->AtomInfo, etc.)*/ } for(a=0;a<I->NAtom;a++) { /* adjust down to point to the count, not the first entry */ if(I->Neighbor[a]>=0) I->Neighbor[a]--; } l=I->Neighbor; } } /*========================================================================*/ CoordSet *ObjectMoleculeChemPyModel2CoordSet(PyObject *model,AtomInfoType **atInfoPtr) { int nAtom,nBond; int a,c; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL,*ai; float *f; BondType *ii,*bond=NULL; int ok=true; int auto_show_lines; int auto_show_nonbonded; int hetatm; int ignore_ids; PyObject *atomList = NULL; PyObject *bondList = NULL; PyObject *atom = NULL; PyObject *bnd = NULL; PyObject *index = NULL; PyObject *crd = NULL; PyObject *tmp = NULL; auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); ignore_ids=!(int)SettingGet(cSetting_preserve_chempy_ids); nAtom=0; nBond=0; if(atInfoPtr) atInfo = *atInfoPtr; atomList = PyObject_GetAttrString(model,"atom"); if(atomList) nAtom = PyList_Size(atomList); else ok=ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't get atom list"); if(ok) { coord=VLAlloc(float,3*nAtom); if(atInfo) VLACheck(atInfo,AtomInfoType,nAtom); } if(ok) { f=coord; for(a=0;a<nAtom;a++) { atom = PyList_GetItem(atomList,a); if(!atom) ok=ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't get atom"); crd = PyObject_GetAttrString(atom,"coord"); if(!crd) ok=ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't get coordinates"); else { for(c=0;c<3;c++) { tmp = PyList_GetItem(crd,c); if (tmp) ok = PConvPyObjectToFloat(tmp,f++); if(!ok) { ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read coordinates"); break; } } } Py_XDECREF(crd); ai = atInfo+a; ai->id = a; /* chempy models are zero-based */ if(!ignore_ids) { if(ok) { /* get chempy atom id if extant */ if(PTruthCallStr(atom,"has","id")) { tmp = PyObject_GetAttrString(atom,"id"); if (tmp) ok = PConvPyObjectToInt(tmp,&ai->id); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read atom identifier"); Py_XDECREF(tmp); } else { ai->id=-1; } } } if(ok) { tmp = PyObject_GetAttrString(atom,"name"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->name,sizeof(AtomName)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read name"); Py_XDECREF(tmp); } if(ok) { if(PTruthCallStr(atom,"has","text_type")) { tmp = PyObject_GetAttrString(atom,"text_type"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->textType,sizeof(TextType)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read text_type"); Py_XDECREF(tmp); } else { ai->textType[0]=0; } } if(ok) { if(PTruthCallStr(atom,"has","vdw")) { tmp = PyObject_GetAttrString(atom,"vdw"); if (tmp) ok = PConvPyObjectToFloat(tmp,&ai->vdw); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read vdw radius"); Py_XDECREF(tmp); } else { ai->vdw=0.0f; } } if(ok) { if(PTruthCallStr(atom,"has","bohr")) { tmp = PyObject_GetAttrString(atom,"bohr"); if (tmp) ok = PConvPyObjectToFloat(tmp,&ai->bohr_radius); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read bohr radius"); Py_XDECREF(tmp); } else { ai->bohr_radius=0.0F; } } if(ok) { if(PTruthCallStr(atom,"has","stereo")) { tmp = PyObject_GetAttrString(atom,"stereo"); if (tmp) ok = PConvPyObjectToInt(tmp,&ai->stereo); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read stereo"); Py_XDECREF(tmp); } else { ai->stereo = 0; } } if(ok) { if(PTruthCallStr(atom,"has","numeric_type")) { tmp = PyObject_GetAttrString(atom,"numeric_type"); if (tmp) ok = PConvPyObjectToInt(tmp,&ai->customType); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read numeric_type"); Py_XDECREF(tmp); } else { ai->customType = cAtomInfoNoType; } } if(ok) { if(PTruthCallStr(atom,"has","formal_charge")) { tmp = PyObject_GetAttrString(atom,"formal_charge"); if (tmp) ok = PConvPyObjectToInt(tmp,&ai->formalCharge); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read formal_charge"); Py_XDECREF(tmp); } else { ai->formalCharge = 0; } } if(ok) { if(PTruthCallStr(atom,"has","partial_charge")) { tmp = PyObject_GetAttrString(atom,"partial_charge"); if (tmp) ok = PConvPyObjectToFloat(tmp,&ai->partialCharge); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read partial_charge"); Py_XDECREF(tmp); } else { ai->partialCharge = 0.0; } } if(ok) { if(PTruthCallStr(atom,"has","flags")) { tmp = PyObject_GetAttrString(atom,"flags"); if (tmp) ok = PConvPyObjectToInt(tmp,(int*)&ai->flags); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read flags"); Py_XDECREF(tmp); } else { ai->flags = 0; } } if(ok) { tmp = PyObject_GetAttrString(atom,"resn"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->resn,sizeof(ResName)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read resn"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"resi"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->resi,sizeof(ResIdent)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read resi"); else ai->resv=AtomResvFromResi(ai->resi); Py_XDECREF(tmp); } if(ok) { if(PTruthCallStr(atom,"has","resi_number")) { tmp = PyObject_GetAttrString(atom,"resi_number"); if (tmp) ok = PConvPyObjectToInt(tmp,&ai->resv); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read resi_number"); Py_XDECREF(tmp); } } if(ok) { tmp = PyObject_GetAttrString(atom,"segi"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->segi,sizeof(SegIdent)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read segi"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"b"); if (tmp) ok = PConvPyObjectToFloat(tmp,&ai->b); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read b value"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"q"); if (tmp) ok = PConvPyObjectToFloat(tmp,&ai->q); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read occupancy"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"chain"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->chain,sizeof(Chain)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read chain"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"hetatm"); if (tmp) ok = PConvPyObjectToInt(tmp,&hetatm); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read hetatm"); else ai->hetatm = hetatm; Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"alt"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->alt,sizeof(Chain)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read alternate conformation"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"symbol"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->elem,sizeof(AtomName)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read symbol"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(atom,"ss"); if (tmp) ok = PConvPyObjectToStrMaxClean(tmp,ai->ssType,sizeof(SSType)-1); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read secondary structure"); Py_XDECREF(tmp); } for(c=0;c<cRepCnt;c++) { atInfo[a].visRep[c] = false; } atInfo[a].visRep[cRepLine] = auto_show_lines; /* show lines by default */ atInfo[a].visRep[cRepNonbonded] = auto_show_nonbonded; /* show lines by default */ if(ok&&atInfo) { AtomInfoAssignParameters(ai); atInfo[a].color=AtomInfoGetColor(ai); } if(!ok) break; } } bondList = PyObject_GetAttrString(model,"bond"); if(bondList) nBond = PyList_Size(bondList); else ok=ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't get bond list"); if(ok) { bond=VLAlloc(BondType,nBond); ii=bond; for(a=0;a<nBond;a++) { bnd = PyList_GetItem(bondList,a); if(!bnd) ok=ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't get bond"); index = PyObject_GetAttrString(bnd,"index"); if(!index) ok=ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't get bond indices"); else { for(c=0;c<2;c++) { tmp = PyList_GetItem(index,c); if (tmp) ok = PConvPyObjectToInt(tmp,&ii->index[c]); if(!ok) { ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read coordinates"); break; } } } if(ok) { tmp = PyObject_GetAttrString(bnd,"order"); if (tmp) ok = PConvPyObjectToInt(tmp,&ii->order); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read bond order"); Py_XDECREF(tmp); } if(ok) { tmp = PyObject_GetAttrString(bnd,"stereo"); if (tmp) ok = PConvPyObjectToInt(tmp,&ii->stereo); else ii->stereo = 0; if(!ok) ii->stereo = 0; Py_XDECREF(tmp); } ii->id=a; if(!ignore_ids) { if(ok) { /* get unique chempy bond id if present */ if(PTruthCallStr(bnd,"has","id")) { tmp = PyObject_GetAttrString(bnd,"id"); if (tmp) ok = PConvPyObjectToInt(tmp,&ii->id); if(!ok) ErrMessage("ObjectMoleculeChemPyModel2CoordSet","can't read bond identifier"); Py_XDECREF(tmp); } else { ii->id=-1; } } } Py_XDECREF(index); ii++; } } Py_XDECREF(atomList); Py_XDECREF(bondList); if(ok) { cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->NTmpBond=nBond; cset->TmpBond=bond; } else { VLAFreeP(bond); VLAFreeP(coord); } if(atInfoPtr) *atInfoPtr = atInfo; if(PyErr_Occurred()) PyErr_Print(); return(cset); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadChemPyModel(ObjectMolecule *I,PyObject *model,int frame,int discrete) { CoordSet *cset = NULL; AtomInfoType *atInfo; int ok=true; int isNew = true; unsigned int nAtom = 0; PyObject *tmp,*mol; if(!I) isNew=true; else isNew=false; if(ok) { if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset=ObjectMoleculeChemPyModel2CoordSet(model,&atInfo); mol = PyObject_GetAttrString(model,"molecule"); if(mol) { if(PyObject_HasAttrString(mol,"title")) { tmp = PyObject_GetAttrString(mol,"title"); if(tmp) { UtilNCopy(cset->Name,PyString_AsString(tmp),sizeof(WordType)); Py_DECREF(tmp); if(!strcmp(cset->Name,"untitled")) /* ignore untitled */ cset->Name[0]=0; } } Py_DECREF(mol); } if(PyObject_HasAttrString(model,"spheroid")&& PyObject_HasAttrString(model,"spheroid_normals")) { tmp = PyObject_GetAttrString(model,"spheroid"); if(tmp) { cset->NSpheroid = PConvPyListToFloatArray(tmp,&cset->Spheroid); if(cset->NSpheroid<0) cset->NSpheroid=0; Py_DECREF(tmp); } tmp = PyObject_GetAttrString(model,"spheroid_normals"); if(tmp) { PConvPyListToFloatArray(tmp,&cset->SpheroidNormal); Py_DECREF(tmp); } } mol = PyObject_GetAttrString(model,"molecule"); nAtom=cset->NIndex; } /* include coordinate set */ if(ok) { if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } else { ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_AllMask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); if(cset->Symmetry&&(!I->Symmetry)) { I->Symmetry=SymmetryCopy(cset->Symmetry); SymmetryAttemptGeneration(I->Symmetry,false,false); } SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadCoords(ObjectMolecule *I,PyObject *coords,int frame) { CoordSet *cset = NULL; int ok=true; int a,l; PyObject *v; float *f; a=0; while(a<I->NCSet) { if(I->CSet[a]) { cset=I->CSet[a]; break; } a++; } if(!PyList_Check(coords)) ErrMessage("LoadsCoords","passed argument is not a list"); else { l = PyList_Size(coords); if (l==cset->NIndex) { cset=CoordSetCopy(cset); f=cset->Coord; for(a=0;a<l;a++) { v=PyList_GetItem(coords,a); /* no error checking */ *(f++)=(float)PyFloat_AsDouble(PyList_GetItem(v,0)); *(f++)=(float)PyFloat_AsDouble(PyList_GetItem(v,1)); *(f++)=(float)PyFloat_AsDouble(PyList_GetItem(v,2)); } } } /* include coordinate set */ if(ok) { if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; SceneCountFrames(); } return(I); } /*========================================================================*/ void ObjectMoleculeBlindSymMovie(ObjectMolecule *I) { CoordSet *frac; int a,c; int x,y,z; float m[16]; if(I->NCSet!=1) { ErrMessage("ObjectMolecule:","SymMovie only works on objects with a single state."); } else if(!I->Symmetry) { ErrMessage("ObjectMolecule:","No symmetry loaded!"); } else if(!I->Symmetry->NSymMat) { ErrMessage("ObjectMolecule:","No symmetry matrices!"); } else if(I->CSet[0]) { frac = CoordSetCopy(I->CSet[0]); CoordSetRealToFrac(frac,I->Symmetry->Crystal); for(x=-1;x<2;x++) for(y=-1;y<2;y++) for(z=-1;z<2;z++) for(a=0;a<I->Symmetry->NSymMat;a++) { if(!((!a)&&(!x)&&(!y)&&(!z))) { c = I->NCSet; VLACheck(I->CSet,CoordSet*,c); I->CSet[c] = CoordSetCopy(frac); CoordSetTransform44f(I->CSet[c],I->Symmetry->SymMatVLA+(a*16)); identity44f(m); m[3] = (float)x; m[7] = (float)y; m[11] = (float)z; CoordSetTransform44f(I->CSet[c],m); CoordSetFracToReal(I->CSet[c],I->Symmetry->Crystal); I->NCSet++; } } frac->fFree(frac); } SceneChanged(); } /*========================================================================*/ void ObjectMoleculeExtendIndices(ObjectMolecule *I) { int a; CoordSet *cs; for(a=-1;a<I->NCSet;a++) { if(a<0) cs=I->CSTmpl; else cs=I->CSet[a]; if(cs) if(cs->fExtendIndices) cs->fExtendIndices(cs,I->NAtom); } } /*========================================================================*/ static CoordSet *ObjectMoleculeMOLStr2CoordSet(char *buffer,AtomInfoType **atInfoPtr) { char *p; int nAtom,nBond; int a,c,cnt,atm,chg; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL; char cc[MAXLINELEN],cc1[MAXLINELEN],resn[MAXLINELEN] = "UNK"; float *f; BondType *ii; BondType *bond=NULL; int ok=true; int auto_show_lines; int auto_show_nonbonded; WordType nameTmp; auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); p=buffer; nAtom=0; if(atInfoPtr) atInfo = *atInfoPtr; /* p=ParseWordCopy(nameTmp,p,sizeof(WordType)-1);*/ p=ParseNCopy(nameTmp,p,sizeof(WordType)-1); p=nextline(p); p=nextline(p); p=nextline(p); if(ok) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&nAtom)!=1) ok=ErrMessage("ReadMOLFile","bad atom count"); } if(ok) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&nBond)!=1) ok=ErrMessage("ReadMOLFile","bad bond count"); } if(ok) { coord=VLAlloc(float,3*nAtom); if(atInfo) VLACheck(atInfo,AtomInfoType,nAtom); } p=nextline(p); /* read coordinates and atom names */ if(ok) { f=coord; for(a=0;a<nAtom;a++) { if(ok) { p=ncopy(cc,p,10); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMOLFile","bad coordinate"); } if(ok) { p=ncopy(cc,p,10); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMOLFile","bad coordinate"); } if(ok) { p=ncopy(cc,p,10); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMOLFile","bad coordinate"); } if(ok) { p=nskip(p,1); p=ncopy(atInfo[a].name,p,3); UtilCleanStr(atInfo[a].name); for(c=0;c<cRepCnt;c++) { atInfo[a].visRep[c] = false; } atInfo[a].visRep[cRepLine] = auto_show_lines; /* show lines by default */ atInfo[a].visRep[cRepNonbonded] = auto_show_nonbonded; /* show lines by default */ } if(ok) { p=nskip(p,2); p=ncopy(cc,p,3); if(sscanf(cc,"%d",&atInfo[a].formalCharge)==1) { if(atInfo[a].formalCharge) { atInfo[a].formalCharge = 4-atInfo[a].formalCharge; } } p=ncopy(cc,p,3); if(sscanf(cc,"%d",&atInfo[a].stereo)!=1) atInfo[a].stereo=0; } if(ok&&atInfo) { atInfo[a].id = a+1; strcpy(atInfo[a].resn,resn); atInfo[a].hetatm=true; AtomInfoAssignParameters(atInfo+a); atInfo[a].color=AtomInfoGetColor(atInfo+a); atInfo[a].alt[0]=0; atInfo[a].segi[0]=0; atInfo[a].resi[0]=0; } p=nextline(p); if(!ok) break; } } if(ok) { bond=VLAlloc(BondType,nBond); ii=bond; for(a=0;a<nBond;a++) { if(ok) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&ii->index[0])!=1) ok=ErrMessage("ReadMOLFile","bad bond atom"); } if(ok) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&ii->index[1])!=1) ok=ErrMessage("ReadMOLFile","bad bond atom"); } if(ok) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&ii->order)!=1) ok=ErrMessage("ReadMOLFile","bad bond order"); } if(ok) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&ii->stereo)!=1) ii->stereo=0; } ii++; if(!ok) break; p=nextline(p); } ii=bond; for(a=0;a<nBond;a++) { ii->index[0]--;/* adjust bond indexs down one */ ii->index[1]--; ii++; } } while(*p) { /* read M CHG records */ p=ncopy(cc,p,6); if(!strcmp(cc,"M CHG")) { p=ncopy(cc,p,3); if(sscanf(cc,"%d",&cnt)==1) { while(cnt--) { p=ncopy(cc,p,4); p=ncopy(cc1,p,4); if(!((*cc)||(*cc1))) break; if((sscanf(cc,"%d",&atm)==1)&& (sscanf(cc1,"%d",&chg)==1)) { atm--; if((atm>=0)&&(atm<nAtom)) atInfo[atm].formalCharge = chg; } } } } p=nextline(p); } if(ok) { cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->NTmpBond=nBond; cset->TmpBond=bond; strcpy(cset->Name,nameTmp); } else { VLAFreeP(bond); VLAFreeP(coord); } if(atInfoPtr) *atInfoPtr = atInfo; return(cset); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadMOLStr(ObjectMolecule *I,char *MOLStr,int frame,int discrete) { int ok = true; CoordSet *cset=NULL; AtomInfoType *atInfo; int isNew; int nAtom; if(!I) isNew=true; else isNew=false; if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset=ObjectMoleculeMOLStr2CoordSet(MOLStr,&atInfo); if(!cset) { ObjectMoleculeFree(I); I=NULL; ok=false; } if(ok) { if(frame<0) frame=I->NCSet; if(I->NCSet<=frame) I->NCSet=frame+1; VLACheck(I->CSet,CoordSet*,frame); nAtom=cset->NIndex; if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } else { ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_MOLMask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadMOLFile(ObjectMolecule *obj,char *fname,int frame,int discrete) { ObjectMolecule* I=NULL; int ok=true; FILE *f; long size; char *buffer,*p; f=fopen(fname,"rb"); if(!f) ok=ErrMessage("ObjectMoleculeLoadMOLFile","Unable to open file!"); else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadMOLFile: Loading from %s.\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); I=ObjectMoleculeReadMOLStr(obj,buffer,frame,discrete); mfree(buffer); } return(I); } /*========================================================================*/ static CoordSet *ObjectMoleculeMOL2Str2CoordSet(char *buffer,AtomInfoType **atInfoPtr,char **next_mol) { char *p; int nAtom,nBond; int a,c; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL; char cc[MAXLINELEN],resn[MAXLINELEN] = "UNK"; float *f; char *last_p; BondType *ii; BondType *bond=NULL; int ok=true; int auto_show_lines; int auto_show_nonbonded; int have_molecule = false; WordType nameTmp; auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); p=buffer; nAtom=0; if(atInfoPtr) atInfo = *atInfoPtr; while((*p)&&ok) { last_p = p; /* top level -- looking for an RTI, assuming p points to the beginning of a line */ p=ParseWordCopy(cc,p,MAXLINELEN); if(cc[0]=='@') { if(WordMatchExact(cc,"@<TRIPOS>MOLECULE",true)||WordMatchExact(cc,"@MOLECULE",true)) { if(have_molecule) { *next_mol = last_p; break; /* next record of multi-mol2 */ } p=ParseNextLine(p); p=ParseNTrim(nameTmp,p,sizeof(WordType)-1); /* get mol name */ p=ParseNextLine(p); p=ParseWordCopy(cc,p,MAXLINELEN); if(sscanf(cc,"%d",&nAtom)!=1) ok=ErrMessage("ReadMOL2File","bad atom count"); else { coord=VLAlloc(float,3*nAtom); if(atInfo) VLACheck(atInfo,AtomInfoType,nAtom); p=ParseNCopy(cc,p,MAXLINELEN); if(sscanf(cc,"%d",&nBond)!=1) { nBond = 0; } } p=ParseNextLine(p); p=ParseNextLine(p); have_molecule=true; } else if(WordMatchExact(cc,"@<TRIPOS>ATOM",true)||WordMatchExact(cc,"@ATOM",true)) { if(!have_molecule) { ok=ErrMessage("ReadMOL2File","@ATOM before @MOLECULE!"); break; } p=ParseNextLine(p); f=coord; for(a=0;a<nAtom;a++) { if(ok) { p=ParseWordCopy(cc,p,MAXLINELEN); if(sscanf(cc,"%d",&atInfo[a].id)!=1) ok=ErrMessage("ReadMOL2File","bad atom id"); } if(ok) { p=ParseWordCopy(cc,p,MAXLINELEN); cc[cAtomNameLen] = 0; if(sscanf(cc,"%s",atInfo[a].name)!=1) ok=ErrMessage("ReadMOL2File","bad atom name"); } if(ok) { p=ParseWordCopy(cc,p,MAXLINELEN); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMOL2File","bad x coordinate"); } if(ok) { p=ParseWordCopy(cc,p,MAXLINELEN); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMOL2File","bad y coordinate"); } if(ok) { p=ParseWordCopy(cc,p,MAXLINELEN); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMOL2File","bad z coordinate"); } if(ok) { p=ParseWordCopy(cc,p,cTextTypeLen); if(sscanf(cc,"%s",atInfo[a].textType)!=1) ok=ErrMessage("ReadMOL2File","bad atom type"); else { /* convert atom type to elem symbol */ char *tt = atInfo[a].textType; char *el = atInfo[a].elem; int elc = 0; while(*tt&&((*tt)!='.')) { *(el++) = *(tt++); elc++; if(elc>cAtomNameLen) break; } *el=0; if(el[2]) el[0] = 0; } } if(ok) { p=ParseWordCopy(cc,p,MAXLINELEN); if(p[0]) { /* skip subst_id */ p=ParseWordCopy(cc,p,MAXLINELEN); if(p[0]) { /* skip subst_name */ p=ParseWordCopy(cc,p,MAXLINELEN); if(p[0]) { if(sscanf(cc,"%f",&atInfo[a].partialCharge)!=1) ok=ErrMessage("ReadMOL2File","bad atom charge"); } } } } p=ParseNextLine(p); for(c=0;c<cRepCnt;c++) { atInfo[a].visRep[c] = false; } atInfo[a].visRep[cRepLine] = auto_show_lines; /* show lines by default */ atInfo[a].visRep[cRepNonbonded] = auto_show_nonbonded; /* show lines by default */ atInfo[a].id = a+1; strcpy(atInfo[a].resn,resn); atInfo[a].hetatm=true; AtomInfoAssignParameters(atInfo+a); atInfo[a].color=AtomInfoGetColor(atInfo+a); atInfo[a].alt[0]=0; atInfo[a].segi[0]=0; atInfo[a].resi[0]=0; } } else if(WordMatchExact(cc,"@<TRIPOS>BOND",true)||WordMatchExact(cc,"@BOND",true)) { if(!have_molecule) { ok=ErrMessage("ReadMOL2File","@ATOM before @MOLECULE!"); break; } p=ParseNextLine(p); if(ok) { bond=VLAlloc(BondType,nBond); ii=bond; for(a=0;a<nBond;a++) { if(ok) { p=ParseWordCopy(cc,p,20); if(sscanf(cc,"%d",&ii->id)!=1) ok=ErrMessage("ReadMOL2File","bad atom id"); } if(ok) { p=ParseWordCopy(cc,p,20); if(sscanf(cc,"%d",&ii->index[0])!=1) ok=ErrMessage("ReadMOL2File","bad source atom id"); } if(ok) { p=ParseWordCopy(cc,p,20); if(sscanf(cc,"%d",&ii->index[1])!=1) ok=ErrMessage("ReadMOL2File","bad target atom id"); } if(ok) { p=ParseWordCopy(cc,p,20); if(!cc[1]) { switch(cc[0]) { case '1': ii->order = 1; break; case '2': ii->order = 2; break; case '3': ii->order = 3; break; } } else if(WordMatchExact("ar",cc,true)) { ii->order = 4; } else if(WordMatchExact("am",cc,true)) { ii->order = 1; } else if(WordMatchExact("un",cc,true)) { ii->order = 1; } else if(WordMatchExact("nc",cc,true)) { ii->order = 0; /* is this legal in PyMOL? */ } else if(WordMatchExact("du",cc,true)) { ii->order = 1; } else ok=ErrMessage("ReadMOL2File","bad bond type"); } ii->stereo = 0; ii++; if(!ok) break; p=ParseNextLine(p); } ii=bond; for(a=0;a<nBond;a++) { ii->index[0]--;/* adjust bond indexs down one */ ii->index[1]--; ii++; } } } else p=ParseNextLine(p); } else p=ParseNextLine(p); } if(ok) { cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->NTmpBond=nBond; cset->TmpBond=bond; strcpy(cset->Name,nameTmp); } else { VLAFreeP(bond); VLAFreeP(coord); } if(atInfoPtr) *atInfoPtr = atInfo; return(cset); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadMOL2Str(ObjectMolecule *I,char *MOLStr,int frame,int discrete) { int ok = true; CoordSet *cset=NULL; AtomInfoType *atInfo; int isNew; int nAtom; char *restart=NULL,*start; int repeatFlag=true; int successCnt = 0; start=MOLStr; while(repeatFlag) { repeatFlag=false; if(!I) isNew=true; else isNew=false; if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } restart=NULL; cset=ObjectMoleculeMOL2Str2CoordSet(start,&atInfo,&restart); if(!cset) { ObjectMoleculeFree(I); I=NULL; ok=false; } if(ok) { if(frame<0) frame=I->NCSet; if(I->NCSet<=frame) I->NCSet=frame+1; VLACheck(I->CSet,CoordSet*,frame); nAtom=cset->NIndex; if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } else { ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_MOLMask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); successCnt++; if(successCnt>1) { if(successCnt==2) { PRINTFB(FB_ObjectMolecule,FB_Actions) " ObjectMolReadMOL2Str: read molecule %d\n",1 ENDFB; } PRINTFB(FB_ObjectMolecule,FB_Actions) " ObjectMolReadMOL2Str: read molecule %d\n",successCnt ENDFB; } } if(restart) { repeatFlag=true; start=restart; frame=frame+1; } } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadMOL2File(ObjectMolecule *obj,char *fname,int frame,int discrete) { ObjectMolecule* I=NULL; int ok=true; FILE *f; long size; char *buffer,*p; f=fopen(fname,"rb"); if(!f) ok=ErrMessage("ObjectMoleculeLoadMOL2File","Unable to open file!"); else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadMOL2File: Loading from %s.\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); I=ObjectMoleculeReadMOL2Str(obj,buffer,frame,discrete); mfree(buffer); } return(I); } /*========================================================================*/ void ObjectMoleculeMerge(ObjectMolecule *I,AtomInfoType *ai, CoordSet *cs,int bondSearchFlag,int aic_mask) { int *index,*outdex,*a2i,*i2a; BondType *bond=NULL; int a,b,c,lb=0,nb,ac,a1,a2; int found; int nAt,nBd,nBond; int expansionFlag = false; AtomInfoType *ai2; int oldNAtom,oldNBond; oldNAtom = I->NAtom; oldNBond = I->NBond; /* first, sort the coodinate set */ index=AtomInfoGetSortedIndex(ai,cs->NIndex,&outdex); for(b=0;b<cs->NIndex;b++) cs->IdxToAtm[b]=outdex[cs->IdxToAtm[b]]; for(b=0;b<cs->NIndex;b++) cs->AtmToIdx[b]=-1; for(b=0;b<cs->NIndex;b++) cs->AtmToIdx[cs->IdxToAtm[b]]=b; ai2=(AtomInfoType*)VLAMalloc(cs->NIndex,sizeof(AtomInfoType),5,true); /* autozero here is important */ for(a=0;a<cs->NIndex;a++) ai2[a]=ai[index[a]]; /* creates a sorted list of atom info records */ VLAFreeP(ai); ai=ai2; /* now, match it up with the current object's atomic information */ for(a=0;a<cs->NIndex;a++) { index[a]=-1; outdex[a]=-1; } c=0; b=0; for(a=0;a<cs->NIndex;a++) { found=false; if(!I->DiscreteFlag) { /* don't even try matching for discrete objects */ lb=b; while(b<I->NAtom) { ac=(AtomInfoCompare(ai+a,I->AtomInfo+b)); if(!ac) { found=true; break; } else if(ac<0) { break; } b++; } } if(found) { index[a]=b; /* store real atom index b for a in index[a] */ b++; } else { index[a]=I->NAtom+c; /* otherwise, this is a new atom */ c++; b=lb; } } /* first, reassign atom info for matched atoms */ /* allocate additional space */ if(c) { expansionFlag=true; nAt=I->NAtom+c; } else { nAt=I->NAtom; } if(expansionFlag) { VLACheck(I->AtomInfo,AtomInfoType,nAt); } /* allocate our new x-ref tables */ if(nAt<I->NAtom) nAt=I->NAtom; a2i = Alloc(int,nAt); i2a = Alloc(int,cs->NIndex); if(nAt) { ErrChkPtr(a2i); } if(cs->NIndex){ ErrChkPtr(i2a); } for(a=0;a<cs->NIndex;a++) /* a is in original file space */ { a1=cs->IdxToAtm[a]; /* a1 is in sorted atom info space */ a2=index[a1]; i2a[a]=a2; /* a2 is in object space */ if(a2<oldNAtom) AtomInfoCombine(I->AtomInfo+a2,ai+a1,aic_mask); else *(I->AtomInfo+a2)=*(ai+a1); } if(I->DiscreteFlag) { if(I->NDiscrete<nAt) { VLACheck(I->DiscreteAtmToIdx,int,nAt); VLACheck(I->DiscreteCSet,CoordSet*,nAt); for(a=I->NDiscrete;a<nAt;a++) { I->DiscreteAtmToIdx[a]=-1; I->DiscreteCSet[a]=NULL; } } I->NDiscrete = nAt; } cs->NAtIndex = nAt; I->NAtom = nAt; FreeP(cs->AtmToIdx); FreeP(cs->IdxToAtm); cs->AtmToIdx = a2i; cs->IdxToAtm = i2a; if(I->DiscreteFlag) { FreeP(cs->AtmToIdx); for(a=0;a<cs->NIndex;a++) { I->DiscreteAtmToIdx[cs->IdxToAtm[a]]=a; I->DiscreteCSet[cs->IdxToAtm[a]] = cs; } } else { for(a=0;a<cs->NAtIndex;a++) cs->AtmToIdx[a]=-1; for(a=0;a<cs->NIndex;a++) cs->AtmToIdx[cs->IdxToAtm[a]]=a; } VLAFreeP(ai); AtomInfoFreeSortedIndexes(index,outdex); /* now find and integrate and any new bonds */ if(expansionFlag) { /* expansion flag means we have introduced at least 1 new atom */ nBond = ObjectMoleculeConnect(I,&bond,I->AtomInfo,cs,bondSearchFlag); if(nBond) { index=Alloc(int,nBond); c=0; b=0; nb=0; for(a=0;a<nBond;a++) { /* iterate over new bonds */ found=false; b=nb; /* pick up where we left off */ while(b<I->NBond) { ac=BondCompare(bond+a,I->Bond+b); if(!ac) { /* zero is a match */ found=true; break; } else if(ac<0) { /* gone past position of this bond */ break; } b++; /* no match yet, keep looking */ } if(found) { index[a]=b; /* existing bond...*/ nb=b+1; } else { /* this is a new bond, save index and increment */ index[a]=I->NBond+c; c++; } } /* first, reassign atom info for matched atoms */ if(c) { /* allocate additional space */ nBd=I->NBond+c; VLACheck(I->Bond,BondType,nBd); for(a=0;a<nBond;a++) /* copy the new bonds */ { a2=index[a]; if(a2 >= I->NBond) { I->Bond[a2] = bond[a]; } } I->NBond=nBd; } FreeP(index); } VLAFreeP(bond); } if(oldNAtom) { if(oldNAtom==I->NAtom) { if(oldNBond!=I->NBond) { ObjectMoleculeInvalidate(I,cRepAll,cRepInvBonds); } } else { ObjectMoleculeInvalidate(I,cRepAll,cRepInvAtoms); } } } /*========================================================================*/ #if 0 ObjectMolecule *ObjectMoleculeLoadPDBFile(ObjectMolecule *obj,char *fname, int frame,int discrete, M4XAnnoType *m4x,PDBInfoRec *pdb_info) { ObjectMolecule *I=NULL; int ok=true; FILE *f; long size; char *buffer,*p; f=fopen(fname,"rb"); if(!f) { PRINTFB(FB_ObjectMolecule,FB_Errors) "ObjectMolecule-ERROR: Unable to open file '%s'\n",fname ENDFB; ok=false; } else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadPDBFile: Loading from %s.\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); I=ObjectMoleculeReadPDBStr(obj,buffer,frame,discrete,m4x); mfree(buffer); } return(I); } #endif /*========================================================================*/ void ObjectMoleculeAppendAtoms(ObjectMolecule *I,AtomInfoType *atInfo,CoordSet *cs) { int a; BondType *ii; BondType *si; AtomInfoType *src,*dest; int nAtom,nBond; if(I->NAtom) { nAtom = I->NAtom+cs->NIndex; VLACheck(I->AtomInfo,AtomInfoType,nAtom); dest = I->AtomInfo+I->NAtom; src = atInfo; for(a=0;a<cs->NIndex;a++) *(dest++)=*(src++); I->NAtom=nAtom; VLAFreeP(atInfo); } else { if(I->AtomInfo) VLAFreeP(I->AtomInfo); I->AtomInfo = atInfo; I->NAtom=cs->NIndex; } nBond=I->NBond+cs->NTmpBond; if(!I->Bond) I->Bond=VLAlloc(BondType,nBond); VLACheck(I->Bond,BondType,nBond); ii=I->Bond+I->NBond; si=cs->TmpBond; for(a=0;a<cs->NTmpBond;a++) { ii->index[0]=cs->IdxToAtm[si->index[0]]; ii->index[1]=cs->IdxToAtm[si->index[1]]; ii->order=si->order; ii->stereo=si->stereo; ii->id=-1; ii++; si++; } I->NBond=nBond; } /*========================================================================*/ CoordSet *ObjectMoleculeGetCoordSet(ObjectMolecule *I,int setIndex) { if((setIndex>=0)&&(setIndex<I->NCSet)) return(I->CSet[setIndex]); else return(NULL); } /*========================================================================*/ void ObjectMoleculeTransformTTTf(ObjectMolecule *I,float *ttt,int frame) { int b; CoordSet *cs; for(b=0;b<I->NCSet;b++) { if((frame<0)||(frame==b)) { cs=I->CSet[b]; if(cs) { if(cs->fInvalidateRep) cs->fInvalidateRep(I->CSet[b],cRepAll,cRepInvCoord); MatrixApplyTTTfn3f(cs->NIndex,cs->Coord,ttt,cs->Coord); } } } } /*========================================================================*/ void ObjectMoleculeSeleOp(ObjectMolecule *I,int sele,ObjectMoleculeOpRec *op) { int a,b,c,s,d,t_i; int a1,ind; float r,rms; float v1[3],v2,*vv1,*vv2,*coord,*vt,*vt1,*vt2; int inv_flag; int hit_flag = false; int ok = true; int cnt; int skip_flag; int match_flag=false; int offset; int priority; CoordSet *cs; AtomInfoType *ai,*ai0,*ai_option; PRINTFD(FB_ObjectMolecule) " ObjectMoleculeSeleOp-DEBUG: sele %d op->code %d\n",sele,op->code ENDFD; if(sele>=0) { SelectorUpdateTable(); /* always run on entry */ switch(op->code) { case OMOP_ALTR: case OMOP_AlterState: PBlock(); /* PBlockAndUnlockAPI() is not safe. * what if "v" is invalidated by another thread? */ break; } /* */ switch(op->code) { case OMOP_AddHydrogens: ObjectMoleculeAddSeleHydrogens(I,sele); break; #ifdef _OLD_CODE if(!ObjectMoleculeVerifyChemistry(I)) { ErrMessage(" AddHydrogens","missing chemical geometry information."); } else { doneFlag=false; while(!doneFlag) { doneFlag=true; a=0; while(a<I->NAtom) { ai=I->AtomInfo + a; s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) if(ObjectMoleculeFillOpenValences(I,a)) { hit_flag=true; doneFlag=false; } a++; /* realize that the atom list may have been resorted */ } } if(hit_flag) { ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); } } break; #endif case OMOP_PrepareFromTemplate: ai0=op->ai; /* template atom */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { ai = I->AtomInfo + a; ai->hetatm=ai0->hetatm; ai->flags=ai0->flags; strcpy(ai->chain,ai0->chain); strcpy(ai->alt,ai0->alt); strcpy(ai->segi,ai0->segi); if(op->i1==1) { /* mode 1, merge residue information */ strcpy(ai->resi,ai0->resi); ai->resv=ai0->resv; strcpy(ai->resn,ai0->resn); } if((ai->elem[0]==ai0->elem[0])&&(ai->elem[1]==ai0->elem[1])) ai->color=ai0->color; else ai->color=AtomInfoGetColor(ai); for(b=0;b<cRepCnt;b++) ai->visRep[b]=ai0->visRep[b]; ai->id=-1; op->i2++; } } break; case OMOP_Pop: for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { if(SelectorMoveMember(s,sele,op->i1)) { op->i2--; op->i3++; } if(!op->i2) break; } } break; case OMOP_PDB1: for(b=0;b<I->NCSet;b++) if(I->CSet[b]) { if((b==op->i1)||(op->i1<0)) for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) ind=I->DiscreteAtmToIdx[a]; else ind=-1; } else ind=I->CSet[b]->AtmToIdx[a]; if(ind>=0) CoordSetAtomToPDBStrVLA(&op->charVLA,&op->i2,I->AtomInfo+a, I->CSet[b]->Coord+(3*ind),op->i3,NULL); op->i3++; } } } break; case OMOP_AVRT: /* average vertex coordinate */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if((priority=SelectorIsMember(s,sele))) { cnt=0; for(b=0;b<I->NCSet;b++) { if(I->CSet[b]) { if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=I->CSet[b]->AtmToIdx[a]; if(a1>=0) { if(!cnt) { VLACheck(op->vv1,float,(op->nvv1*3)+2); VLACheck(op->vc1,int,op->nvv1); } cnt++; vv2=I->CSet[b]->Coord+(3*a1); vv1=op->vv1+(op->nvv1*3); *(vv1++)+=*(vv2++); *(vv1++)+=*(vv2++); *(vv1++)+=*(vv2++); } } } op->vc1[op->nvv1]=cnt; if(cnt) { if(op->vp1) { VLACheck(op->vp1,int,op->nvv1); op->vp1[op->nvv1] = priority; } op->nvv1++; } } } break; case OMOP_StateVRT: /* state vertex coordinate */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if((priority=SelectorIsMember(s,sele))) { cnt=0; b=op->i1; if(b<I->NCSet) if(I->CSet[b]) { if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=I->CSet[b]->AtmToIdx[a]; if(a1>=0) { if(!cnt) { VLACheck(op->vv1,float,(op->nvv1*3)+2); VLACheck(op->vc1,int,op->nvv1); } cnt++; vv2=I->CSet[b]->Coord+(3*a1); vv1=op->vv1+(op->nvv1*3); *(vv1++)+=*(vv2++); *(vv1++)+=*(vv2++); *(vv1++)+=*(vv2++); } } op->vc1[op->nvv1]=cnt; if(cnt) { if(op->vp1) { VLACheck(op->vp1,int,op->nvv1); op->vp1[op->nvv1] = priority; } op->nvv1++; } } } break; case OMOP_SFIT: /* state fitting within a single object */ vt = Alloc(float,3*op->nvv2); cnt = 0; for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { cnt++; break; } } if(cnt) { /* only perform action for selected object */ for(b=0;b<I->NCSet;b++) { rms = -1.0; vt1 = vt; /* reset target vertex pointers */ vt2 = op->vv2; t_i = 0; /* original target vertex index */ if(I->CSet[b]&&(b!=op->i2)) { op->nvv1=0; for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=I->CSet[b]->AtmToIdx[a]; if(a1>=0) { match_flag=false; while(t_i<op->nvv2) { if(op->i1VLA[t_i]==a) {/* same atom? */ match_flag=true; break; } if(op->i1VLA[t_i]<a) { /* catch up? */ t_i++; vt2+=3; } else break; } if(match_flag) { VLACheck(op->vv1,float,(op->nvv1*3)+2); vv2=I->CSet[b]->Coord+(3*a1); vv1=op->vv1+(op->nvv1*3); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); *(vt1++)=*(vt2); *(vt1++)=*(vt2+1); *(vt1++)=*(vt2+2); op->nvv1++; } } } } if(op->nvv1!=op->nvv2) { PRINTFB(FB_Executive,FB_Warnings) "Executive-Warning: Missing atoms in state %d (%d instead of %d).\n", b+1,op->nvv1,op->nvv2 ENDFB; } if(op->nvv1) { if(op->i1!=0) /* fitting flag */ rms = MatrixFitRMS(op->nvv1,op->vv1,vt,NULL,op->ttt); else rms = MatrixGetRMS(op->nvv1,op->vv1,vt,NULL); if(op->i1==2) ObjectMoleculeTransformTTTf(I,op->ttt,b); } else { PRINTFB(FB_Executive,FB_Warnings) "Executive-Warning: No matches found for state %d.\n",b+1 ENDFB; } } VLACheck(op->f1VLA,float,b); op->f1VLA[b]=rms; } VLASize(op->f1VLA,float,I->NCSet); /* NOTE this action is object-specific! */ } FreeP(vt); break; case OMOP_SetGeometry: for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { ai = I->AtomInfo + a; ai->geom=op->i1; ai->valence=op->i2; op->i3++; hit_flag=true; break; } } break; case OMOP_OnOff: for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { hit_flag=true; break; } } break; case OMOP_SaveUndo: /* save undo */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { hit_flag=true; break; } } break; case OMOP_Identify: /* identify atoms */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { VLACheck(op->i1VLA,int,op->i1); op->i1VLA[op->i1++]=I->AtomInfo[a].id; } } break; case OMOP_GetBFactors: ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { op->ff1[op->i1]=ai->b; op->i1++; } ai++; } break; case OMOP_GetOccupancies: ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { op->ff1[op->i1]=ai->q; op->i1++; } ai++; } break; case OMOP_GetPartialCharges: ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { op->ff1[op->i1]=ai->partialCharge; op->i1++; } ai++; } break; case OMOP_IdentifyObjects: /* identify atoms */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { VLACheck(op->i1VLA,int,op->i1); op->i1VLA[op->i1]=I->AtomInfo[a].id; VLACheck(op->obj1VLA,ObjectMolecule*,op->i1); op->obj1VLA[op->i1]=I; op->i1++; } } break; case OMOP_Index: /* identify atoms */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { VLACheck(op->i1VLA,int,op->i1); op->i1VLA[op->i1]=a; /* NOTE: need to incr by 1 before python */ VLACheck(op->obj1VLA,ObjectMolecule*,op->i1); op->obj1VLA[op->i1]=I; op->i1++; } } break; case OMOP_GetObjects: /* identify atoms */ for(a=0;a<I->NAtom;a++) { s=I->AtomInfo[a].selEntry; if(SelectorIsMember(s,sele)) { VLACheck(op->obj1VLA,ObjectMolecule*,op->i1); op->obj1VLA[op->i1]=I; op->i1++; break; } } break; case OMOP_CountAtoms: /* count atoms in object, in selection */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) op->i1++; ai++; } break; case OMOP_PhiPsi: ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { VLACheck(op->i1VLA,int,op->i1); op->i1VLA[op->i1]=a; VLACheck(op->obj1VLA,ObjectMolecule*,op->i1); op->obj1VLA[op->i1]=I; VLACheck(op->f1VLA,float,op->i1); VLACheck(op->f2VLA,float,op->i1); if(ObjectMoleculeGetPhiPsi(I,a,op->f1VLA+op->i1,op->f2VLA+op->i1,op->i2)) op->i1++; } ai++; } break; case OMOP_Cartoon: /* adjust cartoon type */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { ai->cartoon = op->i1; op->i2++; } ai++; } break; case OMOP_Protect: /* protect atoms from movement */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { ai->protekted = op->i1; op->i2++; } ai++; } break; case OMOP_Mask: /* protect atoms from selection */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { ai->masked = op->i1; op->i2++; } ai++; } break; case OMOP_SetB: /* set B-value */ ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { ai->b = op->f1; op->i2++; } ai++; } break; case OMOP_Remove: /* flag atoms for deletion */ ai=I->AtomInfo; if(I->DiscreteFlag) /* for now, can't remove atoms from discrete objects */ ErrMessage("Remove","Can't remove atoms from discrete objects."); else for(a=0;a<I->NAtom;a++) { ai->deleteFlag=false; s=ai->selEntry; if(SelectorIsMember(s,sele)) { ai->deleteFlag=true; op->i1++; } ai++; } break; case OMOP_GetChains: ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { op->ii1[(int)ai->chain[0]]++; op->i1++; } ai++; } break; case OMOP_Spectrum: ai=I->AtomInfo; ai0=NULL; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { skip_flag=false; if(op->i4&&ai0) /* byres and we've done a residue */ if(AtomInfoSameResidue(ai,ai0)) skip_flag=true; if(!skip_flag) { c = (int)(0.49999+op->i1*(op->ff1[op->i3] - op->f1)/op->f2); if(c<0) c=0; if(c>=op->i1) c=op->i1-1; ai->color=op->ii1[c]; /* printf("%8.3 %8.3\n",ai->partial_charge,*/ if(op->i4) { /* byres */ offset = -1; while((a+offset)>=0) { ai0 = I->AtomInfo + a + offset; if(AtomInfoSameResidue(ai,ai0)) { ai0->color = op->ii1[c]; } else break; offset--; } offset = 1; while((a+offset)<I->NAtom) { ai0 = I->AtomInfo + a + offset; if(AtomInfoSameResidue(ai,ai0)) { ai0->color = op->ii1[c]; } else break; offset++; } } ai0=ai; } op->i3++; } ai++; } break; case OMOP_SingleStateVertices: /* same as OMOP_VERT for a single state */ ai=I->AtomInfo; if(op->cs1<I->NCSet) { if(I->CSet[op->cs1]) { b = op->cs1; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { op->i1++; if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=I->CSet[b]->AtmToIdx[a]; if(a1>=0) { VLACheck(op->vv1,float,(op->nvv1*3)+2); vv2=I->CSet[b]->Coord+(3*a1); vv1=op->vv1+(op->nvv1*3); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); op->nvv1++; } } ai++; } } } break; case OMOP_CSetIdxGetAndFlag: ai=I->AtomInfo; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { for(b=op->cs1;b<=op->cs2;b++) { offset = b-op->cs1; if(b<I->NCSet) { if(I->CSet[b]) { if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=I->CSet[b]->AtmToIdx[a]; if(a1>=0) { op->ii1[op->i1*offset+op->i2] = 1; /* presence flag */ vv1=op->vv1+3*(op->i1*offset+op->i2); /* atom-based offset */ vv2=I->CSet[b]->Coord+(3*a1); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); op->nvv1++; } } } } op->i2++; /* atom index field for atoms within selection...*/ } ai++; } break; case OMOP_CSetIdxSetFlagged: ai=I->AtomInfo; hit_flag=false; for(a=0;a<I->NAtom;a++) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { for(b=op->cs1;b<=op->cs2;b++) { offset = b-op->cs1; if(b<I->NCSet) { if(I->CSet[b]) { if(I->DiscreteFlag) { if(I->CSet[b]==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=I->CSet[b]->AtmToIdx[a]; if(a1>=0) { if(op->ii1[op->i1*offset+op->i2]) { /* copy flag */ vv1=op->vv1+3*(op->i1*offset+op->i2); /* atom-based offset */ vv2=I->CSet[b]->Coord+(3*a1); *(vv2++)=*(vv1++); *(vv2++)=*(vv1++); *(vv2++)=*(vv1++); op->nvv1++; hit_flag=true; } } } } } op->i2++; /* atom index field for atoms within selection...*/ } ai++; } break; default: ai = I->AtomInfo; for(a=0;a<I->NAtom;a++) { switch(op->code) { case OMOP_Flag: ai->flags &= op->i2; /* clear flag using mask */ op->i4++; /* no break here is intentional! */ case OMOP_FlagSet: case OMOP_FlagClear: case OMOP_COLR: /* normal atom based loops */ case OMOP_VISI: case OMOP_CheckVis: case OMOP_TTTF: case OMOP_ALTR: case OMOP_LABL: case OMOP_AlterState: s=ai->selEntry; if(SelectorIsMember(s,sele)) { switch(op->code) { case OMOP_Flag: ai->flags |= op->i1; /* set flag */ op->i3++; break; case OMOP_FlagSet: ai->flags |= op->i1; /* set flag */ op->i3++; break; case OMOP_FlagClear: ai->flags &= op->i2; /* clear flag */ op->i3++; break; case OMOP_VISI: if(op->i1<0) for(d=0;d<cRepCnt;d++) ai->visRep[d]=op->i2; else { ai->visRep[op->i1]=op->i2; if(op->i1==cRepCell) I->Obj.RepVis[cRepCell]=op->i2; } break; break; case OMOP_CheckVis: if(ai->visRep[op->i1]) { op->i2 = true; } break; case OMOP_COLR: ai->color=op->i1; op->i2++; break; case OMOP_TTTF: hit_flag=true; break; case OMOP_LABL: if (ok) { if(!op->s1[0]) { ai->label[0]=0; op->i1++; ai->visRep[cRepLabel]=false; hit_flag=true; } else { if(PLabelAtom(&I->AtomInfo[a],op->s1,a)) { op->i1++; ai->visRep[cRepLabel]=true; hit_flag=true; } else ok=false; } } break; case OMOP_ALTR: if (ok) { if(PAlterAtom(&I->AtomInfo[a],op->s1,op->i2,I->Obj.Name,a)) op->i1++; else ok=false; } break; case OMOP_AlterState: if (ok) { if(op->i2<I->NCSet) { cs = I->CSet[op->i2]; if(cs) { if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { if(op->i4) ai_option = I->AtomInfo + a; else ai_option = NULL; if(PAlterAtomState(cs->Coord+(a1*3),op->s1,op->i3,ai_option)) { op->i1++; hit_flag=true; } else ok=false; } } } } break; } break; } break; /* coord-set based properties, iterating only a single coordinate set */ case OMOP_CSetMinMax: case OMOP_CSetCameraMinMax: case OMOP_CSetMaxDistToPt: case OMOP_CSetSumVertices: case OMOP_CSetMoment: if((op->cs1>=0)&&(op->cs1<I->NCSet)) { /* if((I->NCSet==1)&&(SettingGet_i(NULL,I->Obj.Setting,cSetting_static_singletons))) cs=I->CSet[0]; treat static singletons as present in each state else */ cs=I->CSet[op->cs1]; if(cs) { s=ai->selEntry; if(SelectorIsMember(s,sele)) { switch(op->code) { case OMOP_CSetSumVertices: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } add3f(op->v1,coord,op->v1); op->i1++; } break; case OMOP_CSetMinMax: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } if(op->i1) { for(c=0;c<3;c++) { if(*(op->v1+c)>*(coord+c)) *(op->v1+c)=*(coord+c); if(*(op->v2+c)<*(coord+c)) *(op->v2+c)=*(coord+c); } } else { for(c=0;c<3;c++) { *(op->v1+c)=*(coord+c); *(op->v2+c)=*(coord+c); } } op->i1++; } break; case OMOP_CSetCameraMinMax: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } MatrixTransform3f(op->mat1,coord,v1); /* convert to view-space */ coord=v1; if(op->i1) { for(c=0;c<3;c++) { if(*(op->v1+c)>*(coord+c)) *(op->v1+c)=*(coord+c); if(*(op->v2+c)<*(coord+c)) *(op->v2+c)=*(coord+c); } } else { for(c=0;c<3;c++) { *(op->v1+c)=*(coord+c); *(op->v2+c)=*(coord+c); } } op->i1++; } break; case OMOP_CSetMaxDistToPt: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { float dist; coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } dist = (float)diff3f(op->v1,coord); if(dist>op->f1) op->f1=dist; op->i1++; } break; case OMOP_CSetMoment: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { subtract3f(cs->Coord+(3*a1),op->v1,v1); v2=v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]; op->d[0][0] += v2 - v1[0] * v1[0]; op->d[0][1] += - v1[0] * v1[1]; op->d[0][2] += - v1[0] * v1[2]; op->d[1][0] += - v1[1] * v1[0]; op->d[1][1] += v2 - v1[1] * v1[1]; op->d[1][2] += - v1[1] * v1[2]; op->d[2][0] += - v1[2] * v1[0]; op->d[2][1] += - v1[2] * v1[1]; op->d[2][2] += v2 - v1[2] * v1[2]; } break; } } } } break; default: /* coord-set based properties, iterating as all coordsets within atoms */ for(b=0;b<I->NCSet;b++) if(I->CSet[b]) { cs=I->CSet[b]; inv_flag=false; s=ai->selEntry; if(SelectorIsMember(s,sele)) { switch(op->code) { case OMOP_SUMC: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } add3f(op->v1,coord,op->v1); op->i1++; } break; case OMOP_MNMX: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } if(op->i1) { for(c=0;c<3;c++) { if(*(op->v1+c)>*(coord+c)) *(op->v1+c)=*(coord+c); if(*(op->v2+c)<*(coord+c)) *(op->v2+c)=*(coord+c); } } else { for(c=0;c<3;c++) { *(op->v1+c)=*(coord+c); *(op->v2+c)=*(coord+c); } } op->i1++; } break; case OMOP_CameraMinMax: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } MatrixTransform3f(op->mat1,coord,v1); /* convert to view-space */ coord=v1; if(op->i1) { for(c=0;c<3;c++) { if(*(op->v1+c)>*(coord+c)) *(op->v1+c)=*(coord+c); if(*(op->v2+c)<*(coord+c)) *(op->v2+c)=*(coord+c); } } else { for(c=0;c<3;c++) { *(op->v1+c)=*(coord+c); *(op->v2+c)=*(coord+c); } } op->i1++; } break; case OMOP_MaxDistToPt: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { float dist; coord = cs->Coord+3*a1; if(op->i2) /* do we want object-transformed coordinates? */ if(I->Obj.TTTFlag) { transformTTT44f3f(I->Obj.TTT,coord,v1); coord=v1; } dist = (float)diff3f(coord,op->v1); if(dist>op->f1) op->f1=dist; op->i1++; } break; case OMOP_MDST: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { r=(float)diff3f(op->v1,cs->Coord+(3*a1)); if(r>op->f1) op->f1=r; } break; case OMOP_INVA: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) /* selection touches this coordinate set */ inv_flag=true; /* so set the invalidation flag */ break; case OMOP_VERT: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { VLACheck(op->vv1,float,(op->nvv1*3)+2); vv2=cs->Coord+(3*a1); vv1=op->vv1+(op->nvv1*3); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); op->nvv1++; } break; case OMOP_SVRT: /* gives us only vertices for a specific coordinate set */ if(b==op->i1) { if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { VLACheck(op->vv1,float,(op->nvv1*3)+2); VLACheck(op->i1VLA,int,op->nvv1); op->i1VLA[op->nvv1]=a; /* save atom index for later comparisons */ vv2=cs->Coord+(3*a1); vv1=op->vv1+(op->nvv1*3); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); *(vv1++)=*(vv2++); op->nvv1++; } } break; /* Moment of inertia tensor - unweighted - assumes v1 is center of molecule */ case OMOP_MOME: if(I->DiscreteFlag) { if(cs==I->DiscreteCSet[a]) a1=I->DiscreteAtmToIdx[a]; else a1=-1; } else a1=cs->AtmToIdx[a]; if(a1>=0) { subtract3f(cs->Coord+(3*a1),op->v1,v1); v2=v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]; op->d[0][0] += v2 - v1[0] * v1[0]; op->d[0][1] += - v1[0] * v1[1]; op->d[0][2] += - v1[0] * v1[2]; op->d[1][0] += - v1[1] * v1[0]; op->d[1][1] += v2 - v1[1] * v1[1]; op->d[1][2] += - v1[1] * v1[2]; op->d[2][0] += - v1[2] * v1[0]; op->d[2][1] += - v1[2] * v1[1]; op->d[2][2] += v2 - v1[2] * v1[2]; } break; } } switch(op->code) { /* full coord-set based */ case OMOP_INVA: if(inv_flag) { if(op->i1<0) { /* invalidate all representations */ for(d=0;d<cRepCnt;d++) { if(cs->fInvalidateRep) cs->fInvalidateRep(cs,d,op->i2); } } else if(cs->fInvalidateRep) /* invalidate only that particular representation */ cs->fInvalidateRep(cs,op->i1,op->i2); } break; } } /* end coordset section */ break; } ai++; } break; } if(hit_flag) { switch(op->code) { case OMOP_TTTF: ObjectMoleculeTransformTTTf(I,op->ttt,-1); break; case OMOP_LABL: ObjectMoleculeInvalidate(I,cRepLabel,cRepInvText); break; case OMOP_AlterState: /* overly coarse - doing all states, could do just 1 */ if(!op->i3) { /* not read_only? */ ObjectMoleculeInvalidate(I,-1,cRepInvRep); SceneChanged(); } break; case OMOP_CSetIdxSetFlagged: ObjectMoleculeInvalidate(I,-1,cRepInvRep); SceneChanged(); break; case OMOP_SaveUndo: op->i2=true; ObjectMoleculeSaveUndo(I,op->i1,false); break; case OMOP_OnOff: ExecutiveSetObjVisib(I->Obj.Name,op->i1); break; } } /* always run on exit...*/ switch(op->code) { case OMOP_ALTR: case OMOP_AlterState: PUnblock(); break; } /* */ } } /*========================================================================*/ void ObjectMoleculeDescribeElement(ObjectMolecule *I,int index, char *buffer) { AtomInfoType *ai; ai=I->AtomInfo+index; if(ai->alt[0]) sprintf(buffer,"/%s/%s/%s/%s`%s/%s`%s",I->Obj.Name,ai->segi,ai->chain,ai->resn,ai->resi,ai->name,ai->alt); else sprintf(buffer,"/%s/%s/%s/%s`%s/%s",I->Obj.Name,ai->segi,ai->chain,ai->resn,ai->resi,ai->name); } /*========================================================================*/ void ObjectMoleculeGetAtomSele(ObjectMolecule *I,int index, char *buffer) { AtomInfoType *ai; ai=I->AtomInfo+index; if(ai->alt[0]) sprintf(buffer,"/%s/%s/%s/%s`%s/%s`%s",I->Obj.Name,ai->segi,ai->chain,ai->resn,ai->resi, ai->name,ai->alt); else sprintf(buffer,"/%s/%s/%s/%s`%s/%s`",I->Obj.Name,ai->segi,ai->chain,ai->resn,ai->resi, ai->name); } /*========================================================================*/ void ObjectMoleculeGetAtomSeleLog(ObjectMolecule *I,int index, char *buffer,int quote) { AtomInfoType *ai; char quo[5] = ""; if(quote) { quo[0]='"'; quo[1]=0; } if(SettingGet(cSetting_robust_logs)) { ai=I->AtomInfo+index; if(ai->alt[0]) sprintf(buffer,"%s/%s/%s/%s/%s`%s/%s`%s%s",quo,I->Obj.Name,ai->segi,ai->chain,ai->resn,ai->resi, ai->name,ai->alt,quo); else sprintf(buffer,"%s/%s/%s/%s/%s`%s/%s`%s",quo,I->Obj.Name,ai->segi,ai->chain,ai->resn,ai->resi, ai->name,quo); } else { sprintf(buffer,"%s(%s`%d)%s",quo,I->Obj.Name,index+1,quo); } } void ObjectMoleculeGetAtomSeleFast(ObjectMolecule *I,int index, char *buffer) { AtomInfoType *ai; WordType segi,chain,resi,name,alt; ai=I->AtomInfo+index; if(ai->segi[0]) { strcpy(segi,"s;"); strcat(segi,ai->segi); } else { strcpy(segi,"s;''"); } if(ai->chain[0]) { strcpy(chain,"c;"); strcat(chain,ai->chain); } else { strcpy(chain,"c;''"); } if(ai->resi[0]) { strcpy(resi,"i;"); strcat(resi,ai->resi); } else { strcpy(resi,"i;''"); } if(ai->name[0]) { strcpy(name,"n;"); strcat(name,ai->name); } else { strcpy(name,"n;''"); } if(ai->alt[0]) { strcpy(alt,"alt "); strcat(alt,ai->alt); } else { strcpy(alt,"alt ''"); } sprintf(buffer,"(%s&%s&%s&%s&%s&%s)",I->Obj.Name,segi,chain,resi,name,alt); } /*========================================================================*/ int ObjectMoleculeGetNFrames(ObjectMolecule *I) { return I->NCSet; } /*========================================================================*/ void ObjectMoleculeUpdate(ObjectMolecule *I) { int a; OrthoBusyPrime(); for(a=0;a<I->NCSet;a++) if(I->CSet[a]) { OrthoBusySlow(a,I->NCSet); PRINTFD(FB_ObjectMolecule) " ObjectMolecule-DEBUG: updating state %d of \"%s\".\n" , a+1, I->Obj.Name ENDFD; if(I->CSet[a]->fUpdate) I->CSet[a]->fUpdate(I->CSet[a]); } if(I->Obj.RepVis[cRepCell]) { if(I->Symmetry) { if(I->Symmetry->Crystal) { if(I->UnitCellCGO) CGOFree(I->UnitCellCGO); I->UnitCellCGO = CrystalGetUnitCellCGO(I->Symmetry->Crystal); } } } PRINTFD(FB_ObjectMolecule) " ObjectMolecule: updates complete for object %s.\n",I->Obj.Name ENDFD; } /*========================================================================*/ void ObjectMoleculeInvalidate(ObjectMolecule *I,int rep,int level) { int a; PRINTFD(FB_ObjectMolecule) " ObjectMoleculeInvalidate: entered. rep: %d level: %d\n",rep,level ENDFD; if(level>=cRepInvBonds) { VLAFreeP(I->Neighbor); /* set I->Neighbor to NULL */ if(I->Sculpt) { SculptFree(I->Sculpt); I->Sculpt = NULL; } ObjectMoleculeUpdateNonbonded(I); if(level>=cRepInvAtoms) { SelectorUpdateObjectSele(I); } } PRINTFD(FB_ObjectMolecule) " ObjectMoleculeInvalidate: invalidating representations...\n" ENDFD; for(a=0;a<I->NCSet;a++) if(I->CSet[a]) { if(I->CSet[a]->fInvalidateRep) I->CSet[a]->fInvalidateRep(I->CSet[a],rep,level); } PRINTFD(FB_ObjectMolecule) " ObjectMoleculeInvalidate: leaving...\n" ENDFD; } /*========================================================================*/ int ObjectMoleculeMoveAtom(ObjectMolecule *I,int state,int index,float *v,int mode,int log) { int result = 0; CoordSet *cs; if(!(I->AtomInfo[index].protekted==1)) { if(state<0) state=0; if(I->NCSet==1) state=0; state = state % I->NCSet; if((!I->CSet[state])&&(SettingGet_b(I->Obj.Setting,NULL,cSetting_all_states))) state=0; cs = I->CSet[state]; if(cs) { result = CoordSetMoveAtom(I->CSet[state],index,v,mode); cs->fInvalidateRep(cs,cRepAll,cRepInvCoord); } } if(log) { OrthoLineType line,buffer; if(SettingGet(cSetting_logging)) { ObjectMoleculeGetAtomSele(I,index,buffer); sprintf(line,"cmd.translate_atom(\"%s\",%15.9f,%15.9f,%15.9f,%d,%d,%d)\n", buffer,v[0],v[1],v[2],state+1,mode,0); PLog(line,cPLog_no_flush); } } /* if(I->Sculpt) { SculptIterateObject(I->Sculpt,I,state,1); }*/ return(result); } /*========================================================================*/ int ObjectMoleculeInitBondPath(ObjectMolecule *I,ObjectMoleculeBPRec *bp ) { int a; bp->dist = Alloc(int,I->NAtom); bp->list = Alloc(int,I->NAtom); for(a=0;a<I->NAtom;a++) bp->dist[a]=-1; bp->n_atom = 0; return 1; } /*========================================================================*/ int ObjectMoleculePurgeBondPath(ObjectMolecule *I,ObjectMoleculeBPRec *bp ) { FreeP(bp->dist); FreeP(bp->list); return 1; } /*========================================================================*/ int ObjectMoleculeGetBondPaths(ObjectMolecule *I,int atom, int max,ObjectMoleculeBPRec *bp) { /* returns list of bond counts from atom to all others dist and list must be vla array pointers or NULL */ int a,a1,a2,n; int cur; int n_cur; int b_cnt = 0; ObjectMoleculeUpdateNeighbors(I); /* reinitialize dist array (if we've done at least one pass) */ for(a=0;a<bp->n_atom;a++) bp->dist[bp->list[a]]=-1; bp->n_atom = 0; bp->dist[atom] = 0; bp->list[bp->n_atom] = atom; bp->n_atom++; cur = 0; while(1) { b_cnt++; if(b_cnt>max) break; n_cur = bp->n_atom-cur; /* iterate through all current atoms */ if(!n_cur) break; while(n_cur--) { a1 = bp->list[cur++]; n=I->Neighbor[a1]; n++; /* skip cnt */ while(1) { a2=I->Neighbor[n]; n+=2; if(a2<0) break; if(bp->dist[a2]<0) { /* for each atom not yet sampled... */ bp->dist[a2]=b_cnt; bp->list[bp->n_atom]=a2; bp->n_atom++; } } } } return(bp->n_atom); } /*========================================================================*/ int ***ObjectMoleculeGetBondPrint(ObjectMolecule *I,int max_bond,int max_type,int *dim) { int a,b,i,c; int at1,at2; int ***result=NULL; ObjectMoleculeBPRec bp; dim[0]=max_type+1; dim[1]=max_type+1; dim[2]=max_bond+1; result=(int***)UtilArrayMalloc((unsigned int*)dim,3,sizeof(int)); UtilZeroMem(**result,dim[0]*dim[1]*dim[2]*sizeof(int)); ObjectMoleculeInitBondPath(I,&bp); for(a=0;a<I->NAtom;a++) { at1 = I->AtomInfo[a].customType; if((at1>=0)&&(at1<=max_type)) { ObjectMoleculeGetBondPaths(I,a,max_bond,&bp); for(b=0;b<bp.n_atom;b++) { i = bp.list[b]; at2 = I->AtomInfo[i].customType; if((at2>=0)&&(at2<=max_type)) { c=bp.dist[i]; result[at1][at2][c]++; } } } } ObjectMoleculePurgeBondPath(I,&bp); return(result); } /*========================================================================*/ float ObjectMoleculeGetAvgHBondVector(ObjectMolecule *I,int atom,int state,float *v) /* computes average hydrogen bonding vector for an atom */ { float result = 0.0; int a1,a2,n; int vec_cnt = 0; float v_atom[3],v_neigh[3],v_diff[3],v_acc[3] = {0.0,0.0,0.0}; CoordSet *cs; ObjectMoleculeUpdateNeighbors(I); a1 = atom; if(state<0) state=0; if(I->NCSet==1) state=0; state = state % I->NCSet; cs = I->CSet[state]; if(cs) { if(CoordSetGetAtomVertex(cs,a1,v_atom)) { /* atom exists in this C-set */ n=I->Neighbor[atom]; n++; while(1) { a2=I->Neighbor[n]; if(a2<0) break; n+=2; if(I->AtomInfo[a2].protons!=1) { /* ignore hydrogens */ if(CoordSetGetAtomVertex(cs,a2,v_neigh)) { subtract3f(v_atom,v_neigh,v_diff); normalize3f(v_diff); add3f(v_diff,v_acc,v_acc); vec_cnt++; } } } if(vec_cnt) { result = (float)length3f(v_acc); result = result/vec_cnt; normalize23f(v_acc,v); } copy3f(v_acc,v); } } return(result); } /*========================================================================*/ int ObjectMoleculeGetAtomVertex(ObjectMolecule *I,int state,int index,float *v) { int result = 0; if(state<0) state=SettingGet_i(NULL,I->Obj.Setting,cSetting_state)-1; if(state<0) state=SceneGetState(); if(I->NCSet==1) state=0; state = state % I->NCSet; if((!I->CSet[state])&&(SettingGet_b(I->Obj.Setting,NULL,cSetting_all_states))) state=0; if(I->CSet[state]) result = CoordSetGetAtomVertex(I->CSet[state],index,v); return(result); } /*========================================================================*/ int ObjectMoleculeSetAtomVertex(ObjectMolecule *I,int state,int index,float *v) { int result = 0; if(state<0) state=SettingGet_i(NULL,I->Obj.Setting,cSetting_state)-1; if(state<0) state=SceneGetState(); if(I->NCSet==1) state=0; state = state % I->NCSet; if((!I->CSet[state])&&(SettingGet_b(I->Obj.Setting,NULL,cSetting_all_states))) state=0; if(I->CSet[state]) result = CoordSetSetAtomVertex(I->CSet[state],index,v); return(result); } /*========================================================================*/ void ObjectMoleculeRender(ObjectMolecule *I,int state,CRay *ray,Pickable **pick,int pass) { int a; PRINTFD(FB_ObjectMolecule) " ObjectMolecule: rendering %s...\n",I->Obj.Name ENDFD; ObjectPrepareContext(&I->Obj,ray); if(I->UnitCellCGO&&(I->Obj.RepVis[cRepCell])) { if(ray) { CGORenderRay(I->UnitCellCGO,ray,ColorGet(I->Obj.Color), I->Obj.Setting,NULL); } else if(pick&&PMGUI) { } else if(PMGUI) { ObjectUseColor(&I->Obj); CGORenderGL(I->UnitCellCGO,ColorGet(I->Obj.Color), I->Obj.Setting,NULL); } } PRINTFD(FB_ObjectMolecule) " ObjectMolecule: CGO's complete...\n" ENDFD; if(state<0) { for(a=0;a<I->NCSet;a++) if(I->CSet[a]) if(I->CSet[a]->fRender) I->CSet[a]->fRender(I->CSet[a],ray,pick,pass); } else if(state<I->NCSet) { I->CurCSet=state % I->NCSet; if(I->CSet[I->CurCSet]) { if(I->CSet[I->CurCSet]->fRender) I->CSet[I->CurCSet]->fRender(I->CSet[I->CurCSet],ray,pick,pass); } } else if(I->NCSet==1) { /* if only one coordinate set, assume static */ if(SettingGet_b(I->Obj.Setting,NULL,cSetting_static_singletons)) if(I->CSet[0]->fRender) I->CSet[0]->fRender(I->CSet[0],ray,pick,pass); } PRINTFD(FB_ObjectMolecule) " ObjectMolecule: rendering complete for object %s.\n",I->Obj.Name ENDFD; } /*========================================================================*/ void ObjectMoleculeDummyUpdate(ObjectMolecule *I,int mode) { switch(mode) { case cObjectMoleculeDummyOrigin: SceneOriginGet(I->CSet[0]->Coord); break; case cObjectMoleculeDummyCenter: SceneGetPos(I->CSet[0]->Coord); break; } } /*========================================================================*/ ObjectMolecule *ObjectMoleculeDummyNew(int type) { ObjectMolecule *I= NULL; int nAtom; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL; int frame=-1; I=ObjectMoleculeNew(false); nAtom=1; coord=VLAlloc(float,3*nAtom); zero3f(coord); atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->TmpBond=NULL; cset->NTmpBond=0; strcpy(cset->Name,"_origin"); cset->Obj = I; cset->fEnumIndices(cset); ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_IDMask); /* NOTE: will release atInfo */ if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; I->NBond = 0; I->Bond = VLAlloc(BondType,0); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeNew(int discreteFlag) { int a; OOAlloc(ObjectMolecule); ObjectInit((CObject*)I); I->Obj.type=cObjectMolecule; I->NAtom=0; I->NBond=0; I->CSet=VLAMalloc(10,sizeof(CoordSet*),5,true); /* auto-zero */ I->NCSet=0; I->Bond=NULL; I->AtomCounter=-1; I->BondCounter=-1; I->DiscreteFlag=discreteFlag; I->NDiscrete=0; I->UnitCellCGO=NULL; I->Sculpt=NULL; I->CSTmpl=NULL; if(I->DiscreteFlag) { /* discrete objects don't share atoms between states */ I->DiscreteAtmToIdx = VLAMalloc(10,sizeof(int),6,false); I->DiscreteCSet = VLAMalloc(10,sizeof(CoordSet*),5,false); I->NDiscrete=0; } else { I->DiscreteAtmToIdx = NULL; I->DiscreteCSet = NULL; } I->Obj.fRender=(void (*)(struct CObject *, int, CRay *, Pickable **,int))ObjectMoleculeRender; I->Obj.fFree= (void (*)(struct CObject *))ObjectMoleculeFree; I->Obj.fUpdate= (void (*)(struct CObject *)) ObjectMoleculeUpdate; I->Obj.fGetNFrame = (int (*)(struct CObject *)) ObjectMoleculeGetNFrames; I->Obj.fDescribeElement = (void (*)(struct CObject *,int index,char *buffer)) ObjectMoleculeDescribeElement; I->Obj.fGetSettingHandle = (CSetting **(*)(struct CObject *,int state)) ObjectMoleculeGetSettingHandle; I->Obj.fGetCaption = (char *(*)(struct CObject *))ObjectMoleculeGetCaption; I->AtomInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ I->CurCSet=0; I->Symmetry=NULL; I->Neighbor=NULL; for(a=0;a<=cUndoMask;a++) { I->UndoCoord[a]=NULL; I->UndoState[a]=-1; } I->UndoIter=0; return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeCopy(ObjectMolecule *obj) { int a; BondType *i0,*i1; AtomInfoType *a0,*a1; OOAlloc(ObjectMolecule); (*I)=(*obj); I->Symmetry=SymmetryCopy(I->Symmetry); /* null-safe */ I->UnitCellCGO=NULL; I->Neighbor=NULL; I->Sculpt=NULL; for(a=0;a<=cUndoMask;a++) I->UndoCoord[a]=NULL; I->CSet=VLAMalloc(I->NCSet,sizeof(CoordSet*),5,true); /* auto-zero */ for(a=0;a<I->NCSet;a++) { I->CSet[a]=CoordSetCopy(obj->CSet[a]); I->CSet[a]->Obj=I; } if(obj->CSTmpl) I->CSTmpl = CoordSetCopy(obj->CSTmpl); else I->CSTmpl=NULL; I->Bond=VLAlloc(BondType,I->NBond); i0=I->Bond; i1=obj->Bond; for(a=0;a<I->NBond;a++) { *(i0++)=*(i1++); /* copy structure */ } I->AtomInfo=VLAlloc(AtomInfoType,I->NAtom); a0=I->AtomInfo; a1=obj->AtomInfo; for(a=0;a<I->NAtom;a++) *(a0++)=*(a1++); for(a=0;a<I->NAtom;a++) { I->AtomInfo[a].selEntry=0; } return(I); } /*========================================================================*/ void ObjectMoleculeFree(ObjectMolecule *I) { int a; SceneObjectDel((CObject*)I); SelectorPurgeObjectMembers(I); for(a=0;a<I->NCSet;a++) if(I->CSet[a]) { if(I->CSet[a]->fFree) I->CSet[a]->fFree(I->CSet[a]); I->CSet[a]=NULL; } if(I->Symmetry) SymmetryFree(I->Symmetry); VLAFreeP(I->Neighbor); VLAFreeP(I->DiscreteAtmToIdx); VLAFreeP(I->DiscreteCSet); VLAFreeP(I->CSet); VLAFreeP(I->AtomInfo); VLAFreeP(I->Bond); if(I->UnitCellCGO) CGOFree(I->UnitCellCGO); for(a=0;a<=cUndoMask;a++) FreeP(I->UndoCoord[a]); if(I->Sculpt) SculptFree(I->Sculpt); if(I->CSTmpl) if(I->CSTmpl->fFree) I->CSTmpl->fFree(I->CSTmpl); ObjectPurge(&I->Obj); OOFreeP(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadMMDStr(ObjectMolecule *I,char *MMDStr,int frame,int discrete) { int ok = true; CoordSet *cset=NULL; AtomInfoType *atInfo; int isNew; int nAtom; if(!I) isNew=true; else isNew=false; if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset=ObjectMoleculeMMDStr2CoordSet(MMDStr,&atInfo); if(!cset) { VLAFreeP(atInfo); ok=false; } if(ok) { if(!I) I=(ObjectMolecule*)ObjectMoleculeNew(discrete); if(frame<0) frame=I->NCSet; if(I->NCSet<=frame) I->NCSet=frame+1; VLACheck(I->CSet,CoordSet*,frame); nAtom=cset->NIndex; if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; if(cset->fEnumIndices) cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ I->NAtom=nAtom; } else { ObjectMoleculeMerge(I,atInfo,cset,false,cAIC_MMDMask); /* NOTE: will release atInfo */ } if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,false); SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeLoadMMDFile(ObjectMolecule *obj,char *fname, int frame,char *sepPrefix,int discrete) { ObjectMolecule* I=NULL; int ok=true; FILE *f; int oCnt=0; long size; char *buffer,*p; char cc[MAXLINELEN],oName[ObjNameMax]; int nLines; f=fopen(fname,"rb"); if(!f) ok=ErrMessage("ObjectMoleculeLoadMMDFile","Unable to open file!"); else { PRINTFB(FB_ObjectMolecule,FB_Blather) " ObjectMoleculeLoadMMDFile: Loading from %s.\n",fname ENDFB; fseek(f,0,SEEK_END); size=ftell(f); fseek(f,0,SEEK_SET); buffer=(char*)mmalloc(size+255); ErrChkPtr(buffer); p=buffer; fseek(f,0,SEEK_SET); fread(p,size,1,f); p[size]=0; fclose(f); p=buffer; while(ok) { ncopy(cc,p,6); if(sscanf(cc,"%d",&nLines)!=1) break; if(ok) { if(sepPrefix) { I=ObjectMoleculeReadMMDStr(NULL,p,frame,discrete); oCnt++; sprintf(oName,"%s-%02d",sepPrefix,oCnt); ObjectSetName((CObject*)I,oName); ExecutiveManageObject((CObject*)I,true,false); } else { I=ObjectMoleculeReadMMDStr(obj,p,frame,discrete); obj=I; } p=nextline(p); while(nLines--) p=nextline(p); } } mfree(buffer); } return(I); } /*========================================================================*/ ObjectMolecule *ObjectMoleculeReadPDBStr(ObjectMolecule *I,char *PDBStr,int frame, int discrete,M4XAnnoType *m4x,char *pdb_name, char **next_pdb,PDBInfoRec *pdb_info) { CoordSet *cset = NULL; AtomInfoType *atInfo; int ok=true; int isNew = true; unsigned int nAtom = 0; char *start,*restart=NULL; int repeatFlag = true; int successCnt = 0; unsigned int aic_mask = cAIC_PDBMask; SegIdent segi_override=""; /* saved segi for corrupted NMR pdb files */ start=PDBStr; while(repeatFlag) { repeatFlag = false; if(!I) isNew=true; else isNew=false; if(ok) { if(isNew) { I=(ObjectMolecule*)ObjectMoleculeNew(discrete); atInfo = I->AtomInfo; isNew = true; } else { atInfo=VLAMalloc(10,sizeof(AtomInfoType),2,true); /* autozero here is important */ isNew = false; } if(isNew) { AtomInfoPrimeColors(); I->Obj.Color = AtomInfoGetCarbColor(); } cset=ObjectMoleculePDBStr2CoordSet(start,&atInfo,&restart, segi_override,m4x,pdb_name, next_pdb,pdb_info); if(m4x) /* preserve original atom IDs for annotated Metaphorics files */ if(m4x->annotated_flag) aic_mask = (cAIC_b|cAIC_q); nAtom=cset->NIndex; } if(pdb_name&&(*next_pdb)) { /* problematic scenario */ } /* include coordinate set */ if(ok) { if(I->DiscreteFlag&&atInfo) { int a; int fp1 = frame+1; AtomInfoType *ai = atInfo; for(a=0;a<nAtom;a++) { (ai++)->discrete_state = fp1; } } cset->Obj = I; cset->fEnumIndices(cset); if(cset->fInvalidateRep) cset->fInvalidateRep(cset,cRepAll,cRepInvRep); if(isNew) { I->AtomInfo=atInfo; /* IMPORTANT to reassign: this VLA may have moved! */ } else { ObjectMoleculeMerge(I,atInfo,cset,true,aic_mask); /* NOTE: will release atInfo */ } if(isNew) I->NAtom=nAtom; if(frame<0) frame=I->NCSet; VLACheck(I->CSet,CoordSet*,frame); if(I->NCSet<=frame) I->NCSet=frame+1; if(I->CSet[frame]) I->CSet[frame]->fFree(I->CSet[frame]); I->CSet[frame] = cset; if(isNew) I->NBond = ObjectMoleculeConnect(I,&I->Bond,I->AtomInfo,cset,true); if(cset->Symmetry&&(!I->Symmetry)) { I->Symmetry=SymmetryCopy(cset->Symmetry); SymmetryAttemptGeneration(I->Symmetry,false,false); } SceneCountFrames(); ObjectMoleculeExtendIndices(I); ObjectMoleculeSort(I); ObjectMoleculeUpdateIDNumbers(I); ObjectMoleculeUpdateNonbonded(I); successCnt++; if(successCnt>1) { if(successCnt==2){ PRINTFB(FB_ObjectMolecule,FB_Actions) " ObjectMolReadPDBStr: read MODEL %d\n",1 ENDFB; } PRINTFB(FB_ObjectMolecule,FB_Actions) " ObjectMolReadPDBStr: read MODEL %d\n",successCnt ENDFB; } } if(restart) { repeatFlag=true; start=restart; frame=frame+1; } } return(I); } /*========================================================================*/ CoordSet *ObjectMoleculeMMDStr2CoordSet(char *buffer,AtomInfoType **atInfoPtr) { char *p; int nAtom,nBond; int a,c,bPart,bOrder; float *coord = NULL; CoordSet *cset = NULL; AtomInfoType *atInfo = NULL,*ai; char cc[MAXLINELEN]; float *f; BondType *ii,*bond=NULL; int ok=true; int auto_show_lines = (int)SettingGet(cSetting_auto_show_lines); int auto_show_nonbonded = (int)SettingGet(cSetting_auto_show_nonbonded); p=buffer; nAtom=0; if(atInfoPtr) atInfo = *atInfoPtr; if(ok) { p=ncopy(cc,p,6); if(sscanf(cc,"%d",&nAtom)!=1) ok=ErrMessage("ReadMMDFile","bad atom count"); } if(ok) { coord=VLAlloc(float,3*nAtom); if(atInfo) VLACheck(atInfo,AtomInfoType,nAtom); } if(!atInfo) ErrFatal("PDBStr2CoordSet","need atom information record!"); /* failsafe for old version..*/ nBond=0; if(ok) { bond=VLAlloc(BondType,6*nAtom); } p=nextline(p); /* read coordinates and atom names */ if(ok) { f=coord; ii=bond; for(a=0;a<nAtom;a++) { ai=atInfo+a; ai->id=a+1; if(ok) { p=ncopy(cc,p,4); if(sscanf(cc,"%d",&ai->customType)!=1) ok=ErrMessage("ReadMMDFile","bad atom type"); } if(ok) { if(ai->customType<=14) strcpy(ai->elem,"C"); else if(ai->customType<=23) strcpy(ai->elem,"O"); else if(ai->customType<=40) strcpy(ai->elem,"N"); else if(ai->customType<=48) strcpy(ai->elem,"H"); else if(ai->customType<=52) strcpy(ai->elem,"S"); else if(ai->customType<=53) strcpy(ai->elem,"P"); else if(ai->customType<=55) strcpy(ai->elem,"B"); else if(ai->customType<=56) strcpy(ai->elem,"F"); else if(ai->customType<=57) strcpy(ai->elem,"Cl"); else if(ai->customType<=58) strcpy(ai->elem,"Br"); else if(ai->customType<=59) strcpy(ai->elem,"I"); else if(ai->customType<=60) strcpy(ai->elem,"Si"); else if(ai->customType<=61) strcpy(ai->elem,"Du"); else if(ai->customType<=62) strcpy(ai->elem,"Z0"); else if(ai->customType<=63) strcpy(ai->elem,"Lp"); else strcpy(ai->elem,"?"); } for(c=0;c<6;c++) { if(ok) { p=ncopy(cc,p,8); if(sscanf(cc,"%d%d",&bPart,&bOrder)!=2) ok=ErrMessage("ReadMMDFile","bad bond record"); else { if(bPart&&bOrder&&(a<(bPart-1))) { nBond++; ii->index[0]=a; ii->index[1]=bPart-1; ii->order=bOrder; ii->stereo=0; ii->id=-1; ii++; } } } } if(ok) { p=ncopy(cc,p,12); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMMDFile","bad coordinate"); } if(ok) { p=ncopy(cc,p,12); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMMDFile","bad coordinate"); } if(ok) { p=ncopy(cc,p,12); if(sscanf(cc,"%f",f++)!=1) ok=ErrMessage("ReadMMDFile","bad coordinate"); } if(ok) { p=nskip(p,1); p=ncopy(cc,p,5); ai->resv = AtomResvFromResi(cc); sprintf(ai->resi,"%d",ai->resv); } if(ok) { p=nskip(p,6); p=ncopy(cc,p,9); if(sscanf(cc,"%f",&ai->partialCharge)!=1) ok=ErrMessage("ReadMMDFile","bad charge"); } if(ok) { p=nskip(p,10); p=ncopy(cc,p,3); if(sscanf(cc,"%s",ai->resn)!=1) ai->resn[0]=0; ai->hetatm=true; } ai->segi[0]=0; ai->alt[0]=0; if(ok) { p=nskip(p,2); p=ncopy(ai->name,p,4); UtilCleanStr(ai->name); if(ai->name[0]==0) { strcpy(ai->name,ai->elem); sprintf(cc,"%02d",a+1); if((strlen(cc)+strlen(ai->name))>4) strcpy(ai->name,cc); else strcat(ai->name,cc); } for(c=0;c<cRepCnt;c++) { ai->visRep[c] = false; } ai->visRep[cRepLine] = auto_show_lines; /* show lines by default */ ai->visRep[cRepNonbonded] = auto_show_nonbonded; /* show lines by default */ } if(ok) { AtomInfoAssignParameters(ai); ai->color = AtomInfoGetColor(ai); } if(!ok) break; p=nextline(p); } } if(ok) VLASize(bond,BondType,nBond); if(ok) { cset = CoordSetNew(); cset->NIndex=nAtom; cset->Coord=coord; cset->NTmpBond=nBond; cset->TmpBond=bond; } else { VLAFreeP(bond); VLAFreeP(coord); } if(atInfoPtr) *atInfoPtr = atInfo; return(cset); }