c c Sc c -- c c Version 2.1 c c Copyright Michael Lawrence, Biomolecular Research Institute c 343 Royal Parade Parkville Victoria Australia c c This program is designed to compute Sc between two molecules. c It also allows the normal products to be 'merged' into grasp surface c files for display in grasp. c c The algorithm is fully detailed in the paper c c The program includes the freely available subroutines c for calculating Connolly surfaces ( http://www.biohedron.com/msp.html ) c and are provided here according to the conditions of distribution. c c 15/10/99 - integrated into CCP4 suite c The core of the program is unchanged from Mike Lawrence's c original version, however a number of utility s/r's have c been replaced and/or removed. c File opening and reading has been brought into line with c CCP4 standard practice. c c 6/6/00 - fixed up GRASP property handling c stop if there are errors on reading the GRASP input files c c program surface_complementarity c implicit none c logical grasp c include 'defaults.fh' c c CCP4 initialisations call ccpfyp call ccprcs(6, 'SC', '$Date: 2002/06/19 16:19:32 $') c write(6,1) 1 format(' Sc (Version 2.0): A program for determing Shape ', . 'Complementarity',//, . ' Copyright Michael Lawrence, Biomolecular Research ', . 'Institute',/,' 343 Royal Parade Parkville Victoria ', . 'Australia',//,' This program is designed to compute ', . 'Sc between two molecules.',/,' It also allows the ', . 'normal products to be "merged" into grasp surface',/, . ' files for display in grasp.',/) c c write(*,*) '___________________________________________________________________' cman write(*,*) 'INSTRUCTIONS for version modified for CCP4' cman write(*,*) 'Run using:' cman write(*,*) ' sc XYZIN foo.dat [SURFIN1 fooin1 SURFIN2 fooin2' cman write(*,*) ' SURFOUT1 fooout1 SURFOUT2 fooout2]' cman write(*,*) '(where foo.dat is the PDB input file, and the other files' cman write(*,*) 'are optional Grasp surface files) followed by keyworded input.' cman write(*,*) 'The input comprises three sections,' cman write(*,*) 'Section 1. Molecule selection' cman write(*,*) ' Commands for selection are as follows' cman write(*,*) ' MOLECULE N' cman write(*,*) ' CHAIN c - include a particular chain' cman write(*,*) ' ZONE c r c r - include of zone of residues (all atoms thereof)' cman write(*,*) ' AT_EXCL c r a - exclude a particular atom' cman write(*,*) ' AT_INCL c r a - include a particular atom' cman write(*,*) ' where' cman write(*,*) ' N is which molecule to put the subsequent selection in (1 or 2)' cman write(*,*) ' c is the chain id (may be blank, i.e. omitted)' cman write(*,*) ' r is the residue name (not residue type)' cman write(*,*) ' a is the atom name' cman write(*,*) ' Entries for this section appear twice, once for each molecule.' cman write(*,*) 'Section 2. Parameters' cman write(*,*) ' PROBE_RADIUS, used to define the solvent excluded surface ' cman write(*,*) ' DOT_DENSITY, the density of dots on the surface' cman write(*,*) ' TRIM, distance used in generating the peripheral band' cman write(*,*) ' INTERFACE, distance used to determine interface atoms' cman write(*,*) ' WEIGHT, weighting factor in the Sc calculations' cman write(*,*) 'Section 3. Interaction with GRASP' cman write(*,*) ' GRASP_BACKGROUND, value for non-interacting portions of molecule' cman write(*,*) ' GRASP_MATCH , is the tolerance for equivalencing grasp and Sc ' cman write(*,*) ' surface points.' cman write(*,*) ' The Sc value of both surfaces is stored as general ' cman write(*,*) ' property 1 of the Grasp surface file.' cman write(*,*) ' ' cman write(*,*) ' The default values for these parameters (defaults.h) should be ' cman write(*,*) ' suitable for most applications.' cman write(*,*) ' ' cman write(*,*) ' Only atoms within the INTERFACE distance of any "buried" atoms,' cman write(*,*) ' defined in the Connelly routine MDS, are used for calculation of' cman write(*,*) ' the surface complementarity. Subsequently, a periphery band of' cman write(*,*) ' buried points are removed if they lie within a distance TRIM of' cman write(*,*) ' any solvent accesible surface points.' write(*,*) '___________________________________________', . '________________________' write(*,*) ' ' c c Read in atomic radii call sc_radii c c Read in pdb file and assign radii to the protein atoms call readpdb c c Read in atom selection and parameter changes call commands(grasp) c c Write out selection to temporary file call writepdb c c Compute Sc and merge weighted products with grasp surface files (if desired) call sc(grasp) c end c-------------------------------------------------------------------------- subroutine sc_radii c ------------------- c c This subroutine reads in the radius file 'sc_radii' c c The file consists of lines with three fields. The first is the residue type, c the second is the atom name, the third is the radius in angstroms. A wild c card "*" may be used to represent a single character. c c Radius assignment is based on the priority scheme that later radii in the c sc_radii file overwrite earlier wild card assignments. c implicit none c include 'setup.fh' c integer plen,iscrad integer ibeg(200),iend(200),idec(200),ityp(200) real ftoken(200) character*4 ctoken(200) character*80 line c c Initialise unit number for dictionary of radii iscrad = 0 c c Set number of radii counter nrad = 0 c c Open radius file call ccpdpn(iscrad,'SCRADII','READONLY','F',0,0) c c Read file 100 read(iscrad,'(a)',end=200)line c c Use parse to break into fields ntokens = -200 call parse(line,ibeg,iend,ityp,ftoken,ctoken,idec,ntokens) if(ntokens .lt. 3) then write(*,*) 'Error in radius assignment line' write(*,'(a)') line(1:plen(line)) call ccperr(1,' Error in subroutine SC_RADII') end if c c Increment counter and store nrad = nrad + 1 if(nrad.gt.MAXRAD) then write(*,*) 'Too many radii in sc_radii: increase MAXRAD' call ccperr(1,' Error in subroutine SC_RADII') end if c rad_res_type(nrad) = line(ibeg(1):iend(1)) rad_atom_name(nrad) = line(ibeg(2):iend(2)) c if(ityp(3).ne.2) then write(*,*) 'Error in converting string to real' call ccperr(1,' Error in subroutine SC_RADII') end if rad(nrad) = ftoken(3) go to 100 c c File now read. 200 close(iscrad) write(*,250)nrad 250 format(' Number of radii read from sc_radii file:',I5) return c end c----------------------------------------------------------------------------- subroutine readpdb c ------------------ c This subroutine reads in the pdb file c Updated to use CCP4 rwbrook libraries c implicit none include 'setup.fh' character id*4,inscod*1,altcod*1,atnam*4,resnam*3,chnam*1,resno*5 character segid*4 integer i,ifail,ixyzin,iresn,iser,iz,iun integer ichar,jchar,lenstr real occ,biso,u(6) logical found,match external match,xyzinit,xyzopen,xyzadvance,xyzatom,xyzcoord external xyzclose,ccperr,lenstr common /tempfile/iun c c Initialise rwbrook subroutines and open pdb file write(*,fmt='(/,a)') . ' Opening PDB input file with logical name XYZIN' call xyzinit ifail = 0 ixyzin = 0 call xyzopen('XYZIN','INPUT',' ',ixyzin,ifail) c c setup temporary file iun = 0 call ccpdpn(iun,'SCTEMP','SCRATCH','U',0,0) c c read until end of file c n_atoms = 0 100 call xyzadvance(ixyzin,0,0,*101,*101) call xyzatom(ixyzin,iser,atnam,resnam,chnam,iresn,resno, . inscod,altcod,segid,iz,id) n_atoms = n_atoms + 1 if(n_atoms .gt. MAXATOMS) then write(*,*) 'Increase MAXATOMS parameter and re-compile' call ccperr(1,' S/r READPDB: too many atoms in file') end if included(n_atoms) = 0 res_type(n_atoms) = resnam chain(n_atoms) = chnam cpjx Mike points out that the res_name read from the file is right cpjx justified but the comparisions in the ZONE keyword assume left cpjx justification. cpjx Fix is to left-justify the name when storing it cpjx The following construct actually removes _all_ spaces from the cpjx string, but should be sufficient for these purposes res_name(n_atoms) = ' ' ichar = 1 do jchar = 1,lenstr(resno) if (resno(jchar:jchar).ne.' ') then res_name(n_atoms)(ichar:ichar) = resno(jchar:jchar) ichar = ichar + 1 end if end do atom_name(n_atoms) = atnam cpjx not clear what "alt" is - should be reserved? cpjx Later used to check for alternative confirmations (altcod?) alt(n_atoms) = altcod call xyzcoord(ixyzin,'O','U',x(n_atoms),y(n_atoms),z(n_atoms), . occ,biso,u) c c Assign radius found =.false. do i = 1, nrad if(match(res_type(n_atoms),rad_res_type(i)).and. . match(atom_name(n_atoms),rad_atom_name(i))) then rad_atom(n_atoms) = rad(i) found = .true. end if end do if(.not.found) then write(*,*) 'No radius found for' write(*,*) 'Residue ',res_type(n_atoms),' Atom ', . atom_name(n_atoms) call ccperr(1,' S/r READPDB: no radius for residue/atom') end if go to 100 101 write(*,102)n_atoms 102 format(/,' Number of atoms read from PDB file: ',I7) write(*,*) ' ' call xyzclose(ixyzin) return end c c----------------------------------------------------------------------------- c c This subroutine reads the selection commands c c Converted to use CCP4 parser subroutines c subroutine commands(grasp) implicit none logical grasp include 'grasp.fh' include 'setup.fh' logical lmol,match,lwarn character*80 line,command character*1 ch,ch1,ch2 integer mol,nsel,i,i1,i2,plen,iun,nfiles character*80 infile1,infile2,outfile1,outfile2 integer ibeg(200),iend(200),ityp(200),idec(200) real ftoken(200) character*4 ctoken(200) logical lend character*4 a,r,r1,r2 real rp, density, band,sep,weight,dmatch,smin logical ccpexs external ccpexs common /bl16/rp,density,band,sep,weight,dmatch,smin common /tempfile/iun grasp = .false. lmol = .false. lwarn = .false. c c Set the files for Grasp to be the logical names graspin(1) = 'SURFIN1' graspin(2) = 'SURFIN2' graspout(1) = 'SURFOUT1' graspout(2) = 'SURFOUT2' c c Before doing anything else, check to see if any GRASP files c have been assigned c c Check input files nfiles = 0 c call ugtenv('SURFIN1',infile1) if (infile1.ne.' ') then if (ccpexs(infile1)) then nfiles = nfiles + 1 else write(6,*) 'File: ',infile1,' assigned to SURFIN1 not found' call ccperr(2,'SURFIN1 not found') end if end if c call ugtenv('SURFIN2',infile2) if (infile2.ne.' ') then if (ccpexs(infile2)) then nfiles = nfiles + 2 else write(6,*) 'File: ',infile2,' assigned to SURFIN2 not found' call ccperr(2,'SURFIN2 not found') end if end if c c If both input files exist then assume that the user wants c GRASP output if (nfiles.eq.3) then grasp = .true. write(*,*) write(*,*) 'GRASP mode enabled - Grasp output will be produced' write(*,*) write(*,*) 'Grasp input files:' write(*,*) infile1(1:plen(infile1)) write(*,*) infile2(1:plen(infile2)) write(*,*) c Check for output filenames call ugtenv('SURFOUT1',outfile1) if (outfile1.eq.' ') outfile1 = 'SURFOUT1' call ugtenv('SURFOUT2',outfile2) if (outfile2.eq.' ') outfile2 = 'SURFOUT2' write(*,*) 'Grasp output files:' write(*,*) outfile1(1:plen(outfile1)) write(*,*) outfile2(1:plen(outfile2)) write(*,*) else c Grasp interaction is not required/possible write(*,*) write(*,*) + 'GRASP mode disabled - no Grasp output will be produced' write(*,*) end if c c Start of loop to read control input write(*,*)'Selection commands:' write(*,*)'-------------------' 1000 CONTINUE command = ' ' ntokens = 200 line = ' ' C C *********************************************************** CALL PARSER(command, LINE, IBEG, IEND, ITYP, FTOKEN, CTOKEN, IDEC, + NTOKENS, LEND, .false.) C *********************************************************** C if (lend) go to 1001 c if (ntokens.eq.0) go to 1000 c c MOLECULE keyword c ---------------- c if(command .eq. 'MOLE') then if(ntokens .ne. 2) then write(*,*) ' No molecule number given' call ccperr(1,' Syntax: MOLECULE ') else if (ityp(2) .ne. 2) then call ccperr(1,' Error in MOLECULE keyword') end if c nb: no checking here for a proper integer... mol = nint(ftoken(2)) if (mol.ne.1 .and. mol.ne.2) then call ccperr(1,' MOLECULE must be 1 or 2') end if lmol = .true. write(*,10)mol 10 format(/,' Molecule',I2) end if c c CHAIN keyword c ------------- c else if(command.eq.'CHAI'.and.lmol) then if(ntokens .eq. 2) then ch = ctoken(2) nsel = 0 do i = 1, n_atoms if(chain(i) .eq. ch) then nsel = nsel + 1 included(i) = mol end if end do if(nsel .ne. 0) then write(*,20)ch,nsel 20 format(' Number of atoms selected in chain ',A1, . ' = ',I6) else write(*,*)' Chain not found' call ccperr(1,' Error in CHAIN keyword') end if else write(*,*)'Error in CHAIN format' call ccperr(1,' Error in CHAIN keyword') end if c c ZONE keyword c ------------ c else if(command.eq.'ZONE' .and. lmol) then if(ntokens .ne. 3 .and. ntokens .ne. 5 ) then write(*,fmt='(a,/,a)') 'Error in format of ZONE selection', . 'Syntax: ZONE [ ] [ ] ' call ccperr(1,' Error in ZONE keyword') end if if(ntokens .eq. 3) then r1 = ctoken(2) r2 = ctoken(3) ch1 = ' ' ch2 = ' ' else if(ntokens.eq. 5) then r1 = ctoken(3) r2 = ctoken(5) ch1 = ctoken(2) ch2 = ctoken(4) end if c c Find first atom i1 = 0 do i= 1 ,n_atoms - 1 if(match(res_name(i),r1) .and. chain(i) .eq. ch1) then i1 = i go to 50 end if end do 50 if(i1 .eq. 0) then write(*,*) 'First residue in ZONE not found' call ccperr(1,' Error in ZONE keyword') end if c c Finds second atom i2 = 0 do i = i1, n_atoms if(match(res_name(i),r2) .and. chain(i) .eq. ch2) then i2 = i go to 51 end if end do 51 if(i2 .eq. 0) then write(*,*) 'Second residue in ZONE not found' call ccperr(1,' Error in ZONE keyword') end if do i = i1, i2 included(i) = mol end do write(*,35)i2-i1+1 35 format(' Number of atoms selected in zone ',' ',' = ',I6) c c AT_INCL keyword c --------------- c else if(command.eq.'AT_I' .and. lmol) then if(ntokens.eq.3) then ch = ' ' r = ctoken(2) a = ctoken(3) else if(ntokens .eq. 4) then ch = ctoken(2) r = ctoken(3) a = ctoken(4) else write(*,*)' Error in AT_INCL format' call ccperr(1,' Error in AT_INCL keyword') end if i1 = 0 do i = 1, n_atoms if(match(res_name(i),r) .and. match(atom_name(i),a) .and. . chain(i) .eq. ch) then i1 = i included(i) = mol write(*,45)ch,r,a 45 format(' Atom ',A1,' ',A4,' ',A4,' included') go to 52 end if end do 52 if(i1 .eq. 0) then write(*,57)ch,r,a 57 format(' Include atom ',A1,' ',A3,' ',A3,' not found') call ccperr(1,' Error in AT_INCL keyword') end if c c AT_EXCL keyword c --------------- c else if(command .eq. 'AT_E' .and. lmol) then if(ntokens.eq.3) then ch = ' ' r = ctoken(2) a = ctoken(3) else if(ntokens .eq. 4) then ch = ctoken(2) r = ctoken(3) a = ctoken(4) else write(*,*)' Error in AT_EXCL format' call ccperr(1,' Error in AT_EXCL keyword') end if i1 = 0 do i = 1, n_atoms if(match(res_name(i),r) .and. match(atom_name(i),a) .and. . chain(i) .eq. ch) then i1 = i included(i) = 0 write(*,55)ch,r,a 55 format(' Atom ',A1,' ',A4,' ',A4,' excluded') go to 53 end if end do 53 if(i1 .eq. 0) then write(*,58)ch,r,a 58 format(' Exclude atom ',A1,' ',A3,' ',A3,' not found') call ccperr(1,' Error in AT_EXCL keyword') end if c c PROBE_RADIUS keyword c -------------------- c else if(command .eq. 'PROB') then if(ityp(2).ne.2) then write(*,*) 'Error in converting PROBE_RADIUS string to real' call ccperr(1,' Error in PROBE_RADIUS keyword') end if if (ftoken(2).ne.rp) lwarn = .true. rp = ftoken(2) write(*,*) write(*,fmt='(a,f8.2,a)') ' Probe radius set to ',rp, . ' Angstroms' c c TRIM keyword c ------------ c else if(command .eq. 'TRIM') then if(ityp(2).ne.2) then write(*,*) 'Error in converting TRIM string to real' call ccperr(1,' Error in TRIM keyword') end if if (ftoken(2).ne.band) lwarn = .true. band = ftoken(2) write(*,*) write(*,fmt='(a,f8.2,a)') ' Trim width set to ',band, . ' Angstroms' c c INTERFACE keyword c ----------------- c else if (command .eq. 'INTE') then if(ityp(2).ne.2) then write(*,*) 'Error in converting INTERFACE string to real' call ccperr(1,' Error in INTERFACE keyword') end if sep = ftoken(2) write(*,*) write(*,fmt='(a,f8.2,a)')' Distance used to define interface ', . sep,' Angstroms' c c DOT_DENSITY keyword c ------------------- c else if (command .eq. 'DOT_') then if(ityp(2).ne.2) then write(*,*) . 'Error in converting DOT_DENSITY string to real' call ccperr(1,' Error in DOT_DENSITY keyword') end if density = ftoken(2) write(*,*) write(*,fmt='(a,f8.2,a)') ' Dot density for surfaces is ', . density,' per square Angstrom' c c WEIGHT keyword c -------------- c else if (command .eq. 'WEIG') then if(ityp(2).ne.2) then write(*,*) . 'Error in converting WEIGHT string to real' call ccperr(1,' Error in WEIGHT keyword') end if if (ftoken(2).ne.weight) lwarn = .true. weight = ftoken(2) write(*,*) write(*,fmt='(a,f8.2,a)')' Weighting factor is ', . weight,' per square Angstrom' c c GRASP keywords c -------------- c else if (command .eq. 'GRAS') then c command = line(ibeg(1):iend(1)) write(*,fmt='(/,a,a)') ' GRASP command: ',command c if (.not.grasp) then call ccperr(2,'Grasp mode disabled: command ignored') go to 1000 end if c c GRASP_MATCH keyword c ------------------- if (command .eq. 'GRASP_MATCH') then if(ntokens .eq. 2) then if(ityp(2).ne.2) then write(*,*)'Error in converting string to real' call ccperr(1,' Error in GRASP_MATCH keyword') end if dmatch = ftoken(2) write(*,*) ' Grasp maximum vertex matching distance',dmatch else write(*,*) 'Error in GRASP_MATCH format' call ccperr(1,' Error in GRASP_MATCH keyword') end if c c GRASP_BACKGROUND keyword c ------------------------ else if (command .eq. 'GRASP_BACKGROUND') then if(ntokens .eq. 2) then if(ityp(2).ne.2) then write(*,*)'Error in converting string to real' call ccperr(1,' Error in GRASP_BACKGROUND keyword') end if smin = ftoken(2) write(*,*) write(*,*) ' Grasp background Sc value',smin else write(*,*) 'Error in GRASP_BACKGROUND format' call ccperr(1,' Error in GRASP_BACKGROUND keyword') end if end if c c END keyword c ----------- c else if (command(1:3) .eq. 'END') then go to 1001 c c UNRECOGNISED keyword c -------------------- c else write(*,fmt='(/,a,a)')' Unknown command : ',command call ccperr(1,' Unrecognised keyword') end if c go to 1000 c 1001 continue c c Check values of parameters to see if they are still sensible c If values of probe radius or trim have been reset then print a c warning c write(6,fmt='(/,a,/a)')' Parameter values', . ' ----------------' c write(6,fmt='(a,f8.2,a)')' Dot density : ',density, . ' per square Angstrom' if(density .le. 0) then call ccperr (1,' BAD VALUE FOR dot_density: must be gt 0') end if write(6,fmt='(a,f8.2,a)')' Interface separation : ',sep, . ' Angstroms' if(sep .le. 0) then call ccperr (1,' BAD VALUE FOR interface: must be gt 0') end if write(6,fmt='(a,f8.2,a)')' Trim width : ',band, . ' Angstroms' if(band .le. 0) then call ccperr (1,' BAD VALUE FOR trim: must be gt 0') end if write(6,fmt='(a,f8.2,a)')' Probe radius : ',rp, . ' Angstroms' if(rp .le. 0) then call ccperr (1,' BAD VALUE FOR probe_radius: must be gt 0') end if write(6,fmt='(a,f8.2,a)')' Weight factor : ',weight, . ' per square Angstrom' if(weight .le. 0) then call ccperr (1,' BAD VALUE FOR weight: must be gt 0') end if c if(lwarn) then write(6,fmt='(/,5(a,/))') . ' Warning: you have specified one or more of the following', . ' keywords:', . ' PROBE_RADIUS / TRIM / WEIGHT', . ' Cross-comparison of Sc values between proteins is only', . ' valid if the same parameters are used in each calculation!' call ccperr(2, . ' Using non-default values of probe, trim or weight') end if c return end c------------------------------------------------------------------------- subroutine writepdb implicit none include 'setup.fh' integer m1,m2,i,iun common /tempfile/iun write(*,*) write(*,*)'Selected atoms:' write(*,*)'---------------' m1 = 0 m2 = 0 do i= 1, n_atoms if(included(i) .ne. 0) then if(alt(i) .ne. ' ' ) then write(*,fmt='(/,x,a)') . ' Warning: multiple conformations detected:' write(*,fmt='(/,a,a3,a,a4,2(a,a1),a,a4,/)') . ' * Residue ',res_type(i),' Atom ', . atom_name(i),' Altcode ',alt(i),' Chain ', . chain(i),' Resno ',res_name(i) write(*,*)' Multiple conformations not permitted in ', . 'interface atoms -' write(*,*)' these must be removed manually from the PDB', . ' file' call ccperr(1,' Multiple conformations in s/r WRITEPDB') end if if(included(i) .eq. 1) then m1 = m1 + 1 else if(included(i).eq. 2) then m2 = m2 + 1 end if write(iun) x(i),y(i),z(i),rad_atom(i),included(i) end if end do write(*,30) m1 30 format(1x,' Number of atoms for first molecule ',I6) write(*,40) m2 40 format(1x,' Number of atoms for second molecule ',I6) return end c c-------------------------------------------------------------------------- logical function match(c2,c1) c c String c2 is matched to string c1. c c1 may include a trailing * wildcard. If c1 has no wildcard, then c the match must be exact, with both strings having the same number of c non-blank characters. If C1 has a trailing * wildcard, then a match is only c required across equivalent non-blank characters c implicit none character*(*) c1,c2 integer i1,i2,plen i1 = plen(c1) i2 = plen(c2) c c c1 is simply '*', always match if(index(c1,'*').eq.1) then match = .true. return c c No wild cards and strings of unequal length - no match else if(index(c1,'*').eq.0 .and.i1 .ne. i2) then match = .false. return else c c Set length to be equal to last non-wild card character if(index(c1,'*').eq.0) then i1 = plen(c1) else i1 = index(c1,'*') - 1 end if c c C2 must be at least that long if(i2 .lt. i1) then match = .false. return end if c c Now match the characters up until the wildcard if(c1(1:i1).eq.c2(1:i1)) then match= .true. return else match= .false. return endif c end if c end c c---------------------------------------------------------------- integer function plen(str) c -------------------------- c c This function gives the location of the last printable character in c string str character str*(*) integer temp i = len(str) do 100 plen=i,1,-1 temp = IBITS(ichar(str(plen:plen)),0,7) if ((temp.NE.127).AND.(temp.GT.32)) goto 200 100 continue 200 return end c c-------------------------------------------------------------------------- subroutine sc(grasp) C -------------------- c c implicit none include 'big.fh' include 'grasp.fh' logical grasp real atoms(3,MAXATM) real radii(MAXATM) real atmden(MAXATM) integer atten(MAXATM) integer mol(MAXATM) common xp,xn,area,iflag,mi1,mi2 real xp(3,MAXDOTS),area(MAXDOTS),xn(3,MAXDOTS),rn(MAXDOTS) integer*2 iatom,itype,iflag(MAXDOTS) integer mi1(MAXDOTS),mi2(MAXDOTS) integer iun logical doflag(MAXDOTS) integer natom, i, m,ntrim(2),nburied,nblock,j,ndots,m1,m2,it integer nat1,nat2 integer m1_buried,m1_accessible integer m2_buried,m2_accessible real x, y, z, r, dist_min,dis real density,sep,weight,band,rp,dmatch,smin common /bl16/rp,density,band,sep,weight,dmatch,smin common /tempfile/iun integer iundots common /dotfile/iundots real*8 total(2) rewind(iun) c c Read in atom coordinates, radii and mol numbers from scatch file c ----------------------------------------------------------------- natom = 0 200 continue read (iun,end=300) x,y,z,r,m natom = natom + 1 atoms(1,natom) = x atoms(2,natom) = y atoms(3,natom) = z radii(natom) = r mol(natom) = m c c - assign density for this atom atmden(natom) = density go to 200 300 continue close(iun) write (6,320) natom 320 format(1x,' Total number of atoms ', i6) write(*,*) ' ' c c Now determine assign the attention numbers for each atom c -------------------------------------------------------- C nblock = 0 nburied = 0 c - loop all atoms c -------------- nat1 = 0 nat2 = 0 do i = 1,natom if(mol(i) .eq. 1) then nat1 = nat1 + 1 else if(mol(i) .eq. 2) then nat2 = nat2 + 1 end if dist_min=99999. do j = 1,natom c c - find nearest neighbour in other molecule c ---------------------------------------- if(mol(i).ne.mol(j)) then r = dis(atoms(1,i),atoms(1,j)) if(r.lt.dist_min) then dist_min=r end if end if end do c c - check if within separator distance c ---------------------------------- if(dist_min .ge. sep) then c - too far away from other molecule, blocker atom only atten(i) = BLOCKER nblock = nblock + 1 else c - potential interface or neighbouring atom atten(i) = BURIED_FLAGGED nburied = nburied + 1 end if end do c c All assigned write(*,*)'Setting up atoms for surfacing:' write(*,*) ' Potential interface atoms ',nburied write(*,*) ' Blocked atoms ',nblock write(*,*) ' Atoms in Molecule 1 ',nat1 write(*,*) ' Atoms in Molecule 2 ',nat2 write(*,*) write(*,*) 'Output diagnostics from Connolly subroutine MDS:' c c Now compute the surface for the atoms in the interface and its neighbours c ------------------------------------------------------------------------- call mds(rp,natom,atoms,radii,atmden,atten,mol) c c Now read the dot file back in again c ----------------------------------- rewind(iundots) write(*,*)' Surfacing complete.' write(*,*) ndots = 1 m1 = 0 m2 = 0 cpjx added initialisation for m1_buried etc m1_buried = 0 m1_accessible = 0 m2_buried = 0 m2_accessible = 0 666 read(iundots,end=661,err=662) iatom,itype, . xp(1,ndots),xp(2,ndots),xp(3,ndots), . area(ndots), . xn(1,ndots),xn(2,ndots),xn(3,ndots), . iflag(ndots) it = iatom c c Establish two arrays which point to different surfaces c ------------------------------------------------------ if(mol(iatom) .eq. 1) then m1 = m1 + 1 mi1(m1) = ndots if(iflag(ndots).eq.BURIED) then m1_buried = m1_buried + 1 else if(iflag(ndots).eq.ACCESSIBLE) then m1_accessible = m1_accessible + 1 else write(*,*)'problems in dot flavour',iflag(ndots) end if else if(mol(iatom) .eq. 2) then m2 = m2 + 1 mi2(m2) = ndots if(iflag(ndots).eq.BURIED) then m2_buried = m2_buried + 1 else if (iflag(ndots).eq.ACCESSIBLE) then m2_accessible = m2_accessible + 1 else write(*,*)'problems in dot flavour',iflag(ndots) end if else write(*,*) 'problem with molecule number -',mol(iatom) end if ndots = ndots + 1 if (ndots.gt.MAXDOTS) then write(*,fmt='(a,/,a)') . ' Number of surface points in file exceeds MAXDOTS', . ' Increase parameter MAXDOTS and recompile' call ccperr(1, ' Subroutine sc: ndots exceeds MAXDOTS') end if go to 666 661 ndots = ndots-1 close(iundots) write(*,*) 'Statistics of surface point partitioning:' write(*,*) write(*,*) ' Number of dots read = ',ndots write(*,*) ' Number of dots read for molecule 1 = ',m1 write(*,*) ' Number of dots read for molecule 2 = ',m2 write(*,*) ' Number of dots buried for molecule 1 = ', . m1_buried write(*,*) ' Number of dots accessible for molecule 1 = ', . m1_accessible write(*,*) ' Number of dots buried for molecule 2 = ', . m2_buried write(*,*) ' Number of dots accessible for molecule 2 = ', . m2_accessible write(*,*) ' ' c c Cut away the periphery of each surface c -------------------------------------- write(*,*)'Removing peripheral band:' write(*,*) call trim(1,m1,xp,mi1,area,iflag,doflag,total(1),ntrim(1),band) call trim(2,m2,xp,mi2,area,iflag,doflag,total(2),ntrim(2),band) c c Print a message about small buried surface c write(6,fmt='(a,/,a,/,a)') . ' Note: the value of Sc may become meaningless if the buried', . ' area is small, and it will not be a good measure for small', . ' crystal contacts.' c c Compute distance arrays and histograms for each surface c ------------------------------------------------------- c call distance(1,2,xn,xp,mi1,mi2,m1,m2,iflag,doflag,total(1), . ntrim(1),rn,weight) call distance(2,1,xn,xp,mi2,mi1,m2,m1,iflag,doflag,total(2), . ntrim(2),rn,weight) c c Do grasp merging if required c if(grasp) then write(*,*) ' ' write(*,*) ' ' write(*,665) 665 format(' GRASP surface property update:',/) call merge_with_grasp(1,m1,mi1,xp,rn,doflag) call merge_with_grasp(2,m2,mi2,xp,rn,doflag) end if c c Output summary of results call summary c c Finished at last c ---------------- call ccperr (0,'Normal termination') c Error messages c 662 call ccperr(1,' SC: Error in reading dot file') end c c----------------------------------------------------------------------------- c subroutine trim(k,m,xp,mi,area,iflag,doflag,total,ntrim,band) c -------------------------------------------------------------- c implicit none c include 'big.fh' c real xp(3,MAXDOTS),area(MAXDOTS),dis integer*2 iflag(MAXDOTS) integer mi(MAXDOTS) logical doflag(MAXDOTS) real*8 total integer ntrim,i,ip,j,jp,m,k real d,band,nburied,nacc c c Loop over one surface c If a point is buried then see if there is an accessible point within distance band c If so set doflag to false, else retain c nburied = 0 nacc = 0 do i = 1, m ip=mi(i) if(iflag(ip).eq.BURIED) then nburied = nburied + 1 else if(iflag(ip).eq.ACCESSIBLE) then nacc = nacc + 1 end if end do c write(*,1)k,nburied c1 format(' Number of points buried for surface',i2,' is ',1F8.0) c write(*,2)k,nacc c2 format(' Number of points accessible for surface',i2,' is ',1F8.0) c ntrim = 0 total = 0.d0 c do i = 1, m ip = mi(i) doflag(ip) = .false. c - buried point found if(iflag(ip) .eq. BURIED) then doflag(ip) = .true. do j = 1, m jp = mi(j) if(iflag(jp) .eq. ACCESSIBLE) then d = dis(xp(1,ip),xp(1,jp)) if(d .le. band) then doflag(ip) = .false. go to 100 end if end if end do 100 continue c - did this point survive ? c if so accummulate its area if(doflag(ip)) then total = total + area(ip) ntrim = ntrim + 1 end if end if end do c c Print trim statistics c --------------------- c write(*,3) k,ntrim,total 3 format(' Number of points left for molecule',i2, . ' after trim is',i6/ . ' Total area left after trim for this molecule is',f9.2, . ' Angstrom^2.',/) return end c c c----------------------------------------------------------------------- c subroutine distance(k,l,xn,xp,mi,mj,m1,m2,iflag,doflag,total, . ntrim,rn,weight) c --------------------------------------------------------------------- c c This computes the nearest neighbour distance and the accumulates the necssary data c implicit none include 'big.fh' common/SUMM/Dab,Sab real*8 Dab(4),Sab(4) real weight real xp(3,MAXDOTS),xn(3,MAXDOTS),rn(MAXDOTS) integer*2 iflag(MAXDOTS) integer mi(MAXDOTS),mj(MAXDOTS) logical doflag(MAXDOTS) integer*4 neighbour,ntrim,m1,m2 real*8 total,distmin_sum,norm_sum integer i,ip,j,jp,last_bin,istart_bin,k,l,ibin,left_trunc real d,r,distmin,cumperc,c,rleft,rmedian,cumarea,dis,dot,abin,perc integer nbin(0:MAXBIN) integer mbin(-MAXBIN:MAXBIN) logical grasp grasp = .false. c c Zero counters for each histogram bin c ------------------------------------ c c - nbin is the distance histogram c do i = 0, MAXBIN nbin(i) = 0 end do c c - mbin is the normal dot product histogram c do i = -MAXBIN,MAXBIN mbin(i) = 0 end do c c - keep track of the average distance c distmin_sum = 0.d0 c c - keep track of the average normal dot product c norm_sum = 0.d0 c c ------------------------------------------------------------------- c Loop over the entire surface: find and flag neighbour of each point c that we're interested in and store nearest neighbour pointer c ------------------------------------------------------------------- do i = 1, m1 ip = mi(i) c - here we look only at the non-trimmed buried points if(doflag(ip)) then c distmin = 9999999. do j = 1, m2 jp = mj(j) c - for the neighbour we look at all buried points if(iflag(jp) .eq. BURIED) then c - compute distance d = dis(xp(1,ip),xp(1,jp)) c - is this closer than we've seen before if(d .le. distmin) then c - yes, flag neighbour neighbour = jp c - keep nearest distance distmin = d end if end if end do c - having looked at all possible neighbours now accumulate stats distmin_sum = distmin_sum + distmin c - decide which bin to put it into and then add to distance histogram ibin = int(distmin/ binwidth_dist) nbin(ibin) = nbin(ibin) + 1 c - work out dot product jp = neighbour r = dot(xn(1,ip),xn(1,jp)) c - weight dot product cpjx I think the weighting factor is the denominator 2 below? cpjx r = r * exp( - (distmin**2) / 2.) r = r * exp( - (distmin**2) * weight) c - rounding errors a problem, so ensure abs(r) < 1. r = min(0.999,max(r,-0.999)) rn(ip) = r c c darr(ip) = distmin c c - decide on which bin to put it and add to weighted normal histogram ibin = left_trunc(r / binwidth_norm ) mbin(ibin) = mbin(ibin) + 1 c - also accumulate norm_sum = norm_sum + r end if end do c c-------------------------------------------------------------------------------------- c Prepare and dump distance histogram for this surface c ---------------------------------------------------- write(*,1) k,l 1 format(/, . '$TABLE : Distance between surfaces D(',i1,'->',i1,') :',/, . '$GRAPHS :Area:N:1,3::Cumul_Area:N:1,4::Percentage:N:1,5:', . ':Cumul_Percentage:N:1,6:$$',/, . ' From To Area Cumulative_Area Perentage ', . ' Cumulative_Percentage $$ $$',/) c c Determine the last distance bin that has anything in it c last_bin = 0 do i = MAXBIN, 0, -1 if(nbin(i) .ne. 0) then last_bin = i go to 30 end if end do 30 continue c c Accumulate percentages and area from all filled distance bins c cumperc = 0. cumarea = 0. c c Loop all active bins c -------------------- do i = 0, last_bin c - abin has the total area associated with points in the i'th distance bin c - (for this we assume that each point has approximately the same area) abin = total * nbin(i) / ntrim cumarea = cumarea + abin c - perc is the percentage of the total area in this bin perc = abin * 100. / total c - accumulate the total percentage area thus far c = cumperc + perc c - check for the 50th percentile and interpolate between the bin limits if( cumperc .le. 50. .and . c .ge. 50.) then rleft = i * binwidth_dist c - rmedian then gives us the median distance rmedian = rleft + (50. - cumperc) * binwidth_dist / . (c - cumperc) end if c - store accumulated percentage cumperc = c c - dump detail for this bin write(*,2) i*binwidth_dist,(i+1)*binwidth_dist,abin,cumarea, . perc,cumperc 2 format(f5.2,' ',f5.2,t14,f7.1,t23,f10.1,t41,f7.1,t59,f7.1) end do write(*,3) 3 format('$$') c c Dump means and medians on distance c ---------------------------------- write(*,22) k,l,distmin_sum/ntrim,k,l,rmedian 22 format(/' Mean separation: surface ',i1,' to surface ',i1,' = ', . f6.3,' Angstroms',/, . ' Median separation: surface ',i1,' to surface ',i1,' = ', . f6.3,' Angstroms') Dab(k)=distmin_sum/ntrim Dab(k+2)=rmedian c c-------------------------------------------------------------------------------------- c Prepare and dump weighted normal histogram for this surface c ----------------------------------------------------------- c write(*,6) k,l,k,l 6 format(// . '$TABLE : Surface complementarity S(',i1,'->',i1,') :',/, . '$GRAPHS :S(',i1,'->',i1,'):N:1,3: :Percentage:N:1,4: $$',/, . ' From To Number Percentage ', . ' $$ $$',/) c c Determine last bin used for weighted normal product c --------------------------------------------------- last_bin = 0 do i = MAXBIN, 0, -1 if(mbin(i) .ne. 0) then last_bin = i go to 31 end if end do c c Determine first bin used for weighted normal product c --------------------------------------------------- istart_bin = 0 31 do i = -MAXBIN,0 if(mbin(i) .ne. 0) then istart_bin = i go to 32 end if end do 32 continue c c Accumulate percentage for medians c --------------------------------- cumperc = 0. c c Loop all active bins c -------------------- do i = istart_bin, last_bin c - accumulate percentage perc = mbin(i) * 100. / ntrim c = cumperc + perc c - check for 50th percentile, interpolate for median if( cumperc .le. 50. .and . c .ge. 50.) then rleft = i * binwidth_norm rmedian = rleft + (50. - cumperc) * binwidth_norm / . (c - cumperc) end if cumperc = c c - dump bin info write(*,61) -i*binwidth_norm,-(i+1)*binwidth_norm, mbin(i),perc 61 format(f6.2,' ',f6.2,t17,i7,t32,f5.1) end do write(*,4) 4 format('$$') c c - Dump overall weighted normal product histogram stats c ---------------------------------------------------- c write(*,23) k,l,-norm_sum/ntrim,k,l,-rmedian 23 format(/' Mean normal product:' . ,' surface ',i1,' to surface ',i1,' = ',f6.3/ . ' Median normal product:', . ' surface ',i1,' to surface ',i1,' = ',f6.3) Sab(k)=-norm_sum/ntrim Sab(k+2)=-rmedian c c Zero norm sums c c do i = 1, 3 c sum_norm(i) = 0.d0 c end do c c do i = 1, np(k) c if( doflag(i,k) ) then c sum_norm(1) = sum_norm(1) + xn(i,k) c sum_norm(2) = sum_norm(2) + yn(i,k) c sum_norm(3) = sum_norm(3) + zn(i,k) c end if c end do c c dsum_norm = 0. c do j = 1, 3 c dsum_norm = dsum_norm + sum_norm(j) ** 2 c end do c c dsum_norm = sqrt(dsum_norm) c c do j = 1, 3 c mean_norm(j) = sum_norm(j) / dsum_norm c end do c c sum_dot = 0.d0 c c do i = 1, np(k) c if ( doflag(i,k) ) then c sum_dot = sum_dot + xn(i,k) * mean_norm(1) c . + yn(i,k) * mean_norm(2) c . + zn(i,k) * mean_norm(3) c end if c end do c c sum_dot = sum_dot / ntrim(k) c c write(*,1) k, sum_dot c 1 format(/' For surface',i2,' the mean normal spread is ',f6.3/) c return end c----------------------------------------------------------------- c integer function left_trunc(r) c ------------------------------ c if(r .ge. 0) then left_trunc = int(r) else left_trunc = int(r) - 1 end if return end c c-------------------------------------------------------------------------- c c c molecular dot surface computation subroutine: c c Nb This subroutine copyright Michael Connolly c obtain the original version from http://www.biohedron.com/msp.html c c The following changes have been made: c 1. "stop" statements have been replaced with calls to ccperr c 2. s/r "cross" is assumed to be the subroutine cross from the ccp4 c modlib library. This has a different ordering of arguments to that c used originally in mds. c 3. "open" statement has been replaced with a call to ccpdpn, and c unit number is in the variable "iundots", which is exported to c the calling routine via the "dotfile" common block. c No other changes have been made. c subroutine mds(rp,natom,atoms,radii,atmden,atten,mol) real rp integer natom integer iundots real atoms(3,natom) real radii(natom) real atmden(natom) integer atten(natom) integer mol(natom) c MAXATM parameter should have the same value in the calling program integer MAXATM parameter (MAXATM=15000) c MAXNBR parameter should have the same value in the putpnt subroutine integer MAXNBR, MAXSUB, MAXPRB, MAXLOW parameter (MAXNBR=200) parameter (MAXSUB=100000) parameter (MAXPRB=150000) parameter (MAXLOW=120000) integer i, k, l, m, n, mm, mc integer i1, i2, i3, i4, nnbr, nbur (MAXATM), nbb, jnbr integer knbr, lnbr integer ilat, nlat, nlow, nnear, nprobe, iprb, iprb2 integer isign, is0, isub, nsub, ipnt, npnt real contain real pradius, bb, bb2, radmax, d, d2, density, bridge real ri, rj, rk, eri, erj, erk, erl real rij, rci, rcj, rb, area, rad real dij, djk, dik, asymm, dmin, dt, wijk, swijk, far real anglei, anglej, dtq, e, rs, cs, ps, ts real hijk, dtijk2, rkp2, dm, edens real dot real triple real dis2 real dis real disptl c neighbor arrays: integer nbr(MAXNBR) real disnbr(MAXNBR),nbrstd(MAXNBR) logical nbrusd(MAXNBR) c probe placement arrays: real heights(MAXPRB), probes(3,MAXPRB), alts(3,MAXPRB) integer pi1(MAXPRB), pi2(MAXPRB), pi3(MAXPRB) integer lowprbs(MAXLOW),nears(MAXLOW) c subdivision arrays: real subs(3,MAXSUB),pnts(3,MAXSUB) real lats(3,MAXSUB) c accessibility array: logical access(MAXATM) integer is(3) real uij(3),tij(3),uik(3),tik(3),uijk(3) real utb(3),bijk(3),pijk(3),pij(3) real pi(3),pj(3),axis(3), qij(3),qji(3) real north(3) real pqi(3),pqj(3),o(3),south(3) real eqvec(3),cen(3),pcen(3),vp(3,3),vector(3,3) real vtemp(3),vql(3) logical pcusp, between c common block arrays: integer nsp(3) real burco(3,MAXNBR,MAXATM),burrad(MAXNBR,MAXATM),pprev(3) logical bprev common /mdsblk/nsp,pradius,nbur,burco,burrad,pprev,bprev common /dotfile/iundots c check for atom overflow if (natom .gt. MAXATM) call ccperr(1, 'mds: too many atoms') c initialization pradius = rp do 25 k = 1,3 nsp(k) = 0 25 continue nprobe = 0 do 50 i = 1,natom access(i) = .false. 50 continue do 60 k = 1,3 pprev(k) = 1000000.0 o(k) = 0.0 60 continue radmax = 0.0 do 70 i = 1,natom if (radii(i) .gt. radmax) radmax = radii(i) 70 continue bb = 4 * radmax + 4 * rp bb2 = bb * bb c open output file iundots = 0 call ccpdpn(iundots,'SCDOTS','SCRATCH','U',0,0) cpjx open (unit=7,file='dots', cpjx :form='unformatted',status='scratch',err=9020) c main loop: do 1000 i1 = 1, natom if (atten(i1) .le. 0) go to 1000 ri = radii(i1) eri = ri + rp density = atmden(i1) nnbr = 0 nbur(i1) = 0 nbb = 0 do 100 i2 = 1,natom if (i1 .eq. i2) go to 100 if (atten(i2) .le. 0) go to 100 if (mol(i2) .eq. mol(i1)) then d2 = dis2(atoms(1,i1),atoms(1,i2)) if (d2 .le. 0.0) go to 9040 bridge = ri + radii(i2) + 2 * rp if (d2 .ge. bridge * bridge) go to 100 nnbr = nnbr + 1 if (nnbr .gt. MAXNBR) go to 9060 nbr(nnbr) = i2 else if (atten(i1) .lt. 5) go to 100 d2 = dis2(atoms(1,i1),atoms(1,i2)) if (d2 .lt. bb2) nbb = nbb + 1 bridge = ri + radii(i2) + 2 * rp if (d2 .ge. bridge * bridge) go to 100 nbur(i1) = nbur(i1) + 1 if (nbur(i1) .gt. MAXNBR) go to 9060 do 90 k = 1,3 burco (k, nbur(i1),i1) = atoms(k,i2) 90 continue burrad(nbur(i1),i1) = radii(i2) end if 100 continue if (atten(i1) .eq. 6 .and. nbb .le. 0) go to 1000 if (nnbr .le. 0) then access(i1) = .true. go to 620 end if c sort neighbors by distance from atom i1 do 110 n = 1,nnbr i2 = nbr(n) disnbr(n) = dis2(atoms(1,i1),atoms(1,i2)) nbrusd(n) = .false. 110 continue do 120 n = 1,nnbr l = 0 dmin = 1000000.0 do 115 m = 1,nnbr if (nbrusd(m)) go to 115 if (disnbr(m) .lt. dmin) then dmin = disnbr(m) l = m end if 115 continue if (l .le. 0) call ccperr(1, 'insertion sort fails') nbrstd(n) = nbr(l) nbrusd(l) = .true. 120 continue c second loop: do 600 jnbr = 1,nnbr i2 = nbr(jnbr) if (i2 .le. i1) go to 600 rj = radii(i2) erj = rj + rp density = (atmden(i1) + atmden(i2))/2 dij = dis(atoms(1,i1),atoms(1,i2)) do 140 k = 1,3 uij(k) = (atoms(k,i2) - atoms(k,i1))/dij 140 continue asymm = (eri * eri - erj * erj) / dij between = (abs(asymm) .lt. dij) do 160 k = 1,3 tij(k) = 0.5 * (atoms(k,i1) + atoms(k,i2)) + 0.5 * asymm * uij(k) 160 continue far = (eri + erj) * (eri + erj) - dij * dij if (far .le. 0.0) go to 600 far = sqrt (far) contain = dij * dij - (ri - rj) * (ri - rj) if (contain .le. 0.0) go to 600 contain = sqrt (contain) rij = 0.5 * far * contain / dij if (nnbr .le. 1) then access(i1) = .true. access(i2) = .true. go to 320 end if c third loop: do 300 knbr = 1,nnbr i3 = nbr(knbr) if (i3 .le. i2) go to 300 rk = radii (i3) erk = rk + rp djk = dis(atoms(1,i2),atoms(1,i3)) if (djk .ge. erj + erk) go to 300 dik = dis(atoms(1,i1), atoms(1,i3)) if (dik .ge. eri + erk) go to 300 if (atten(i1) .le. 1 .and. atten(i2) .le. 1 .and. :atten(i3) .le. 1) go to 300 do 220 k = 1,3 uik(k) = (atoms(k,i3) - atoms(k,i1)) / dik 220 continue dt = dot(uij, uik) if (dt .ge. 1.0) go to 255 if (dt .le. -1.0) go to 255 wijk = acos(dt) if (wijk .le. 0.0) go to 255 swijk = sin(wijk) if (swijk .le. 0.0) go to 255 call cross (uijk,uij,uik) do 230 k = 1,3 uijk(k) = uijk(k) / swijk 230 continue call cross (utb,uijk,uij) asymm = (eri * eri - erk * erk) / dik do 235 k = 1,3 tik(k) = 0.5 * (atoms(k,i1) + atoms(k,i3)) + 0.5 * asymm * uik(k) 235 continue dt = 0.0 do 240 k = 1,3 dt = dt + uik(k) * (tik(k) - tij(k)) 240 continue do 250 k = 1,3 bijk(k) = tij(k) + utb(k) * dt / swijk 250 continue hijk = eri * eri - dis2 (bijk, atoms (1,i1)) if (hijk .gt. 0.0) go to 257 c collinear and other 255 continue dtijk2 = dis2(tij,atoms(1,i3)) rkp2 = erk * erk - rij * rij if (dtijk2 .lt. rkp2) go to 600 go to 300 257 continue hijk = sqrt(hijk) do 280 is0 = 1,2 isign = 3 - 2 * is0 do 260 k = 1,3 pijk(k) = bijk(k) + isign * hijk * uijk(k) 260 continue c check for collision do 270 lnbr = 1,nnbr i4 = nbrstd(lnbr) erl = radii(i4) + rp if (i4 .eq. i2) go to 270 if (i4 .eq. i3) go to 270 if (dis2(pijk,atoms(1,i4)) .le. erl * erl) go to 280 270 continue c new probe position nprobe = nprobe + 1 if (nprobe .gt. MAXPRB) go to 9100 if (isign .gt. 0) then pi1(nprobe) = i1 pi2(nprobe) = i2 pi3(nprobe) = i3 else pi1(nprobe) = i2 pi2(nprobe) = i1 pi3(nprobe) = i3 end if heights(nprobe) = hijk do 275 k = 1,3 probes(k, nprobe) = pijk(k) alts(k,nprobe) = isign * uijk(k) 275 continue access(i1) = .true. access(i2) = .true. access(i3) = .true. 280 continue 300 continue 320 continue c return to second loop c toroidal surface generation if (atten(i1) .le. 1 .and. atten(i2) .le. 1) go to 600 if (rp .le. 0.0) go to 600 density = (atmden(i1) + atmden(i2))/2 rci = rij * ri / eri rcj = rij * rj / erj rb = rij - rp if (rb .le. 0.0) rb = 0.0 rs = (rci + 2 * rb + rcj) / 4 e = rs/rij edens = e * e * density call subcir (tij, rij, uij, edens, MAXSUB, subs, nsub, ts) if (nsub .le. 0) go to 600 do 400 isub = 1, nsub do 325 k = 1,3 pij(k) = subs(k, isub) 325 continue c check for collision do 330 lnbr = 1, nnbr i4 = nbrstd (lnbr) if (i4 .eq. i2) go to 330 erl = radii(i4) + rp if (dis2 (pij, atoms(1,i4)) .le. erl * erl) go to 400 330 continue c no collision, toroidal arc generation access(i1) = .true. access(i2) = .true. if (atten(i1) .eq. 6 .and. atten(i2) .eq. 6 .and. nbur(i1) .le. 0) :go to 400 do 340 k = 1, 3 pi(k) = (atoms(k, i1) - pij(k)) / eri pj(k) = (atoms(k, i2) - pij(k)) / erj 340 continue call cross (axis, pi, pj) call normal (axis) dtq = rp * rp - rij * rij pcusp = (dtq .gt. 0.0 .and. between) if (pcusp) then c point cusp -- two shortened arcs dtq = sqrt (dtq) do 350 k = 1, 3 qij(k) = tij(k) - dtq * uij(k) qji(k) = tij(k) + dtq * uij(k) 350 continue do 360 k = 1, 3 pqi(k) = (qij(k) - pij(k))/rp pqj(k) = (qji(k) - pij(k))/rp 360 continue else c no cusp do 370 k = 1, 3 pqi(k) = pi(k) + pj(k) 370 continue call normal (pqi) do 375 k = 1, 3 pqj(k) = pqi(k) 375 continue end if dt = dot (pi, pqi) if (dt .ge. 1.0) go to 400 if (dt .le. -1.0) go to 400 anglei = acos (dt) dt = dot (pqj, pj) if (dt .ge. 1.0) go to 400 if (dt .le. -1.0) go to 400 anglej = acos (dt) c convert two arcs to points if (atten(i1) .ge. 2) then call subarc (pij, rp, axis, density, pi, pqi, MAXSUB, pnts, npnt, . ps) do 380 ipnt = 1, npnt area = ps * ts * disptl (tij, uij, pnts (1, ipnt)) / rij call putpnt (2, i1, atten (i1), pnts (1, ipnt), area, pij, . atoms (1, i1)) 380 continue end if if (atten(i2) .ge. 2) then call subarc (pij, rp, axis, density, pqj, pj, MAXSUB, pnts, npnt, . ps) do 390 ipnt = 1, npnt area = ps * ts * disptl (tij, uij, pnts (1, ipnt)) / rij call putpnt (2, i2, atten (i2), pnts (1, ipnt) ,area, pij, . atoms(1,i2)) 390 continue end if 400 continue 600 continue c return to first loop 620 continue if (.not. access(i1)) go to 1000 if (atten(i1) .le. 1) go to 1000 if (atten(i1) .eq. 6 .and. nbur(i1) .le. 0) go to 1000 density = atmden (i1) c convex surface generation if (nnbr .le. 0) then north (1) = 0.0 north(2) = 0.0 north(3) = 1.0 south(1) = 0.0 south(2) = 0.0 south(3) = -1.0 eqvec(1) = 1.0 eqvec(2) = 0.0 eqvec(3) = 0.0 else i2 = nbrstd(1) do 630 k = 1,3 north(k) = atoms(k, i1) - atoms(k, i2) 630 continue call normal(north) vtemp(1) = north(2) * north(2) + north(3) * north(3) vtemp(2) = north(1) * north(1) + north(3) * north(3) vtemp(3) = north(1) * north(1) + north(2) * north(2) call normal(vtemp) dt = dot (vtemp, north) if (abs(dt) .gt. 0.99) then vtemp(1) = 1.0 vtemp(2) = 0.0 vtemp(3) = 0.0 end if call cross (eqvec, north, vtemp) call normal (eqvec) call cross (vql, eqvec, north) rj = radii (i2) erj = rj + rp dij = dis (atoms(1, i1), atoms (1, i2)) do 635 k = 1,3 uij(k) = (atoms(k, i2) - atoms (k, i1)) / dij 635 continue asymm = (eri * eri - erj * erj) / dij do 640 k = 1, 3 tij(k) = 0.5 * (atoms (k, i1) + atoms (k, i2)) + 0.5 * asymm . * uij (k) 640 continue far = (eri + erj) * (eri + erj) - dij * dij if (far .le. 0.0) call ccperr(1, 'imaginary far') far = sqrt (far) contain = dij * dij - (ri - rj) * (ri - rj) if (contain .le. 0.0) call ccperr(1, 'imaginary contain') contain = sqrt (contain) rij = 0.5 * far * contain / dij do 645 k = 1,3 pij(k) = tij(k) + rij * vql(k) south(k) = (pij(k) - atoms(k, i1)) / eri 645 continue if (triple(north, south, eqvec) .le. 0.0) . call ccperr(1, 'non-positive frame') end if call subarc (o, ri, eqvec, density, north, south, MAXSUB, lats, . nlat, cs) c debugging print statement: c write (6, 656) i1, nlat, nnbr, nbur, nbb 656 format (1x, 'i1, nlat, nnbr, nbur, nbb:', 5i6) if (nlat .le. 0) go to 1000 do 800 ilat = 1, nlat c project onto north vector dt = dot (lats (1, ilat), north) do 660 k = 1, 3 cen(k) = atoms(k, i1) + dt * north(k) 660 continue rad = ri * ri - dt * dt if (rad .le. 0.0) go to 800 rad = sqrt (rad) call subcir (cen, rad, north, density, MAXSUB, pnts, npnt, ps) if (npnt .le. 0) go to 800 area = cs * ps do 750 ipnt = 1, npnt do 680 k = 1, 3 pcen(k) = atoms(k, i1) + (eri/ri) * (pnts(k,ipnt) - atoms(k,i1)) 680 continue c check for collision do 690 lnbr = 2, nnbr i4 = nbrstd (lnbr) erl = radii (i4) + rp if (dis2 (pcen, atoms(1, i4)) .le. erl * erl) go to 750 690 continue c no collision, put point call putpnt (1, i1, atten(i1), pnts(1,ipnt), area, pcen, . atoms(1,i1)) 750 continue 800 continue 1000 continue c end of first loop c write (6, 1010) nprobe c1010 format (1x, 'nprobe = ', i8) if (rp .le. 0.0) go to 1600 c concave surface generation c collect low probes nlow = 0 do 1050 iprb = 1, nprobe if (heights(iprb) .ge. rp) go to 1050 nlow = nlow + 1 if (nlow .gt. MAXLOW) call ccperr(1, 'low probe overflow') lowprbs (nlow) = iprb 1050 continue c write (6, 1060) nlow c1060 format (1x, 'number of low probes = ', i4) c probe loop: do 1500 iprb = 1, nprobe is(1) = pi1 (iprb) is(2) = pi2 (iprb) is(3) = pi3 (iprb) do 1120 k = 1, 3 pijk(k) = probes (k, iprb) uijk(k) = alts (k, iprb) 1120 continue hijk = heights (iprb) i1 = is(1) i2 = is(2) i3 = is(3) if (atten(i1) .eq. 6 .and. atten(i2) .eq. 6 .and. :atten(i3) .eq. 6 .and. nbur(i1) .le. 0) go to 1500 density = (atmden(i1) + atmden(i2) + atmden(i3)) / 3 c gather nearby low probes nnear = 0 if (nlow .le. 0) go to 1135 do 1130 l = 1, nlow iprb2 = lowprbs(l) if (iprb2 .eq. iprb) go to 1130 d2 = dis2 (pijk, probes (1, iprb2)) if (d2 .ge. 4 * rp * rp) go to 1130 nnear = nnear + 1 if (nnear .gt. MAXLOW) call ccperr(1, 'near low overflow') nears (nnear) = iprb2 1130 continue 1135 continue c set up vectors from probe center to three atoms do 1160 m = 1, 3 i = is (m) do 1140 k = 1,3 vp(k,m) = atoms(k, i) - pijk(k) 1140 continue call normal (vp(1,m)) 1160 continue c set up vectors normal to three cutting planes call cross (vector (1,1), vp(1,1), vp (1,2)) call cross (vector (1,2), vp(1,2), vp (1,3)) call cross (vector (1,3), vp(1,3), vp (1,1)) call normal (vector(1,1)) call normal (vector(1,2)) call normal (vector(1,3)) c find latitude of highest vertex of triangle dm = -1.0 mm = 0 do 1180 m = 1, 3 dt = dot (uijk, vp (1,m)) if (dt .gt. dm) then dm = dt mm = m end if 1180 continue c create arc for selecting latitudes do 1200 k = 1,3 south(k) = - uijk(k) 1200 continue call cross (axis, vp (1, mm), south) call normal (axis) call subarc (o, rp, axis, density, vp (1,mm), south, MAXSUB, lats, . nlat, cs) if (nlat .le. 0) go to 1500 c probe latitude loop: do 1400 ilat = 1, nlat dt = dot (lats(1,ilat),south) do 1250 k = 1,3 cen(k) = dt * south(k) 1250 continue rad = rp * rp - dt * dt if (rad .le. 0.0) go to 1400 rad = sqrt (rad) call subcir (cen, rad, south, density, MAXSUB, pnts, npnt, ps) if (npnt .le. 0) go to 1400 area = cs * ps do 1350 ipnt = 1, npnt c check against 3 planes do 1270 m = 1,3 dt = dot (pnts (1, ipnt), vector (1,m)) if (dt .ge. 0.0) go to 1350 1270 continue do 1280 k = 1,3 pnts(k, ipnt) = pijk(k) + pnts(k, ipnt) 1280 continue c if low probe, check for inside nearby low probe if (hijk .ge. rp) go to 1320 if (nnear .le. 0) go to 1320 do 1300 n = 1, nnear iprb2 = nears(n) d2 = dis2 (pnts (1, ipnt), probes (1, iprb2)) if (d2 .lt. rp * rp) go to 1350 1300 continue 1320 continue c determine which atom the surface point is closest to mc = 0 dmin = 2 * rp do 1330 m = 1,3 i = is(m) d = dis (pnts(1,ipnt), atoms (1,i)) - radii(i) if (d .lt. dmin) then dmin = d mc = m end if 1330 continue i = is (mc) call putpnt (3, i, atten(i), pnts (1, ipnt), area, pijk, . atoms (1, i)) 1350 continue 1400 continue 1500 continue c end of probe loop c branch for van der Waals surface: 1600 continue c finished writing surface points c close (7) write (6, 2050) nsp(1), nsp(2), nsp(3), nsp(1)+nsp(2)+nsp(3) 2050 format (1x, ' Number of surface points: ',/, . ' Convex ',i8,/ . ' Toroidal ',i8,/ . ' Concave ',i8,/ . ' Total ',i8,/) c c stop if no atoms found in interface if ((nsp(1)+nsp(2)+nsp(3)).eq.0) then call ccperr(1, ' No atoms found in interface.') endif return c error branches: cpjx 9020 call ccperr(1, 'mds: cannot open output file') 9040 call ccperr(1, 'mds: coincident atoms') 9060 call ccperr(1, 'mds: neighbor overflow') 9100 call ccperr(1, 'mds: probe overflow') end c subroutines c put point to disk: subroutine putpnt (itype, i, atten, coor, area, pcen, acen) integer itype, i, atten real area real coor(3), pcen(3), acen(3) integer k, lbur, nbur real erl, pradius real outnml(3) c MAXNBR should have the same value here and in the mds subroutine integer MAXNBR parameter (MAXNBR=200) parameter (MAXATM=15000) integer nsp(3) real burco(3,MAXNBR,MAXATM),burrad(MAXNBR,MAXATM),pprev(3) logical bprev common /mdsblk/nsp, pradius, nbur(MAXATM), burco, burrad, pprev, . bprev logical buried integer*2 i2, itype2, iflag integer iundots common /dotfile/iundots real dis2 buried = .false. i2 = i itype2 = itype go to (100, 200, 300, 400, 400, 400) atten call ccperr(1, 'mds: invalid atom attention number in putpnt') 100 continue return 200 continue write(iundots) i2, itype2, coor 250 format (1x,i5,i2,3f8.3) go to 900 300 continue write(iundots) i2, itype2, coor, area 350 format (1x, i5, i2, 4f8.3) go to 900 400 continue c calculate outward pointing unit normal vector if (pradius .le. 0.0) then do 410 k = 1, 3 outnml(k) = coor(k) - acen(k) 410 continue call normal(outnml) else do 420 k = 1,3 outnml(k) = (pcen(k) - coor(k))/pradius 420 continue end if if (atten .ge. 5) go to 500 write(iundots) i2, itype2, coor, area, outnml 450 format (1x, i5, i2, 7f8.3) go to 900 500 continue c determine whether buried c first check whether probe changed c if (dis2 (pcen, pprev) .le. 0.0) then buried = bprev else c check for collisions with neighbors in other molecules do 530 lbur = 1, nbur(i2) erl = burrad (lbur,i2) + pradius if (dis2 (pcen, burco (1, lbur,i2)) .le. erl * erl) then buried = .true. go to 540 end if 530 continue 540 continue do 545 k = 1, 3 pprev(k) = pcen(k) 545 continue end if bprev = buried if (atten .eq. 6) go to 600 if (buried) then iflag = 1 else iflag = 0 end if write(iundots) i2, itype2, coor, area, outnml, iflag 550 format (1x, i5, i2, 7f12.8, i2) go to 900 600 continue if (.not. buried) return write(iundots) i2, itype2, coor, area, outnml 900 continue nsp(itype) = nsp(itype) + 1 return end c elementary mathematical subroutines: real function triple (c, d, e) real c(3), d(3), e(3) real cd(3) real dot call cross (cd, c, d) triple = dot (cd, e) return end real function dis2 (c, d) real c(3), d(3) integer k real s, cd s = 0.0 do 100 k = 1, 3 cd = c(k) - d(k) s = s + cd * cd 100 continue dis2 = s return end real function dis (c, d) real c(3), d(3) real s real dis2 real sqrt s = dis2 (c, d) if (s .lt. 0.0) s = 0.0 dis = sqrt (s) return end subroutine normal (x) real x(3) integer k real s s = 0.0 do 100 k = 1, 3 s = s + x(k) * x(k) 100 continue if (s .le. 0.0) call ccperr(1, 'zero vec in normal') s = sqrt (s) do 200 k = 1, 3 x(k) = x(k) / s 200 continue return end c distance from point to line: real function disptl (cen, axis, pnt) real cen(3), axis(3), pnt(3) real vec(3) integer k real dt, d2 real dot real sqrt do 100 k = 1, 3 vec(k) = pnt(k) - cen(k) 100 continue dt = dot (vec, axis) d2 = vec(1) * vec(1) + vec(2) * vec(2) + vec(3) * vec(3) - dt * dt if (d2 .le. 0.0) d2 = 0.0 disptl = sqrt (d2) return end c subdivision routines: subroutine subdiv (cen, rad, x, y, angle, density, maxsub, pnts, . npnt, ps) integer maxsub, npnt real rad, angle, density, ps real cen(3) real x(3), y(3) real pnts (3, maxsub) real a, c, s, delta integer n, k, i real sqrt delta = 1.0 / (sqrt(density) * rad) a = - delta / 2 n = 0 do 100 i = 1, maxsub a = a + delta if (a .gt. angle) go to 150 n = n + 1 c = rad * cos (a) s = rad * sin (a) do 50 k = 1, 3 pnts(k, i) = cen(k) + c * x(k) + s * y(k) 50 continue 100 continue if (a + delta .lt. angle) call ccperr(1, 'too many subdivisions') 150 continue npnt = n if (npnt .gt. 0) then ps = rad * angle / npnt else ps = 0.0 end if return end subroutine subarc (cen, rad, axis, density, x, v, maxsub, pnts, . npnt, ps) integer maxsub, npnt real rad, density, ps real cen(3), axis(3), x(3), v(3) real pnts (3, maxsub) real y(3) real angle real dt1, dt2 real dot real atan2 call cross (y, axis, x) dt1 = dot (v, x) dt2 = dot (v, y) angle = atan2 (dt2, dt1) if (angle .lt. 0.0) angle = angle + 2 * 3.1415926535 call subdiv (cen, rad, x, y, angle, density, maxsub, pnts, npnt, . ps) return end subroutine subcir (cen, rad, axis, density, maxsub, pnts, npnt, . ps) integer maxsub, npnt real rad, density, ps real cen(3), axis(3) real pnts(3, maxsub) real v1(3), v2(3), x(3), y(3) real angle, dt real dot c arbitrary unit vector perpendicular to axis v1(1) = axis(2) * axis(2) + axis(3) * axis(3) v1(2) = axis(1) * axis(1) + axis(3) * axis(3) v1(3) = axis(1) * axis(1) + axis(2) * axis(2) call normal (v1) dt = dot(v1,axis) if (abs (dt) .gt. 0.99) then v1(1) = 1.0 v1(2) = 0.0 v1(3) = 0.0 end if call cross (v2, axis, v1) call normal (v2) call cross (x, axis, v2) call normal (x) call cross (y, axis, x) angle = 2.0 * 3.1415926535 call subdiv (cen, rad, x, y, angle, density, maxsub, pnts, npnt, . ps) return end c c Copyright 1986 by Michael L. Connolly c All rights reserved c------------------------------------------------------------------------------- c subroutine merge_with_grasp(ngrasp,m,mi,xp,rn,doflag) c implicit none include 'grasp.fh' real rp, density, band,sep,weight,dmatch,smin common /bl16/rp,density,band,sep,weight,dmatch,smin real xp(3,*),rn(*) logical doflag(*) integer mi(*),ngrasp,m c real buffer(MAXVERT) integer ibuffer4(MAXVERT),i,nvert,icurv,iprop,plen,iform,istart integer igraspin,igraspout integer icomma,ngrid,ntri,l,ip,iblank equivalence (buffer,ibuffer4) external plen real vert(3,MAXVERT),curv(MAXVERT),dmin,d,s integer*4 ibuff integer nmcitm,result integer*2 ibuffer(3) equivalence (ibuffer,buffer) character*80 line,prop(10),pline character*1 c1 logical done character*20 gp(MAXPROP) logical have_prop(MAXPROP) integer igprop c c Set up grasp property order gp(1) = 'potential' gp(2) = 'curvature' gp(3) = 'distance' gp(4) = 'gproperty' gp(5) = 'g2property' gp(6) = 'vertexcolor' igprop = 4 c c We have none to start with do i = 1, MAXPROP have_prop(i) = .false. end do c Test grasp file for compatibility c This is necessary because files written on SGs might not be c readable on other systems with different "endedness" c c At Mike Lawrence's suggestion (PJX implemented): c read in the first integer*4 of the surface file using qread c This should have value = 80 c c If not then assume that we have an incompatibly written file c and abort the merging procedure igraspin = 0 write(*,fmt='(/,x,a,i2,a,/)') + 'Testing GRASP file',ngrasp,' for compatibility...' call qopen(igraspin,graspin(ngrasp),'READONLY') call qmode(igraspin,6,nmcitm) result = 0 call qread(igraspin,ibuff,1,result) call qclose(igraspin) if (ibuff.eq.80 .and. result.eq.0) then write(*,*) write(*,*) 'Grasp file okay...proceeding' else write(*,fmt='(/,x,a,a8,/,x,a,/)') + 'Problems reading Grasp file ',graspin(ngrasp), + 'Maybe it was created on a different system?' call ccperr(2, ' Aborting merge with Grasp surface') return end if c Proceed with merging Grasp surface and sc values igraspin = 0 igraspout = 0 call ccpdpn(igraspin,graspin(ngrasp),'READONLY','U',0,0) call ccpdpn(igraspout,graspout(ngrasp),'UNKNOWN','U',0,0) write(*,*) write(*,40)ngrasp 40 format( ' Header information from grasp file',I2) write(*,*) '------------------------------------' c write(*,*) read(igraspin) line write(igraspout) line write(*,*) line(1:plen(line)) if(index(line,'format=1').ne.0) then iform=1 else if(index(line,'format=2').ne.0) then iform=2 else call ccperr(1, 'Unknown format for GRASP input file') end if read(igraspin) line write(igraspout) line write(*,*) 'Keywords=',line(1:plen(line)) read(igraspin) pline write(*,*) 'Properties=',pline(1:plen(pline)) write(*,*) write(*,*) 'Decoding property list....' C Now determine which properties are present in the file and handle accordingly iprop = 0 iblank = index(pline,' ') cpjx Bizarre bug - using the rwbrook routines seems to cause different cpjx values to be returned for iblank. The check on the length of cpjx pline seems to work around this... if(plen(pline) .eq. 0 .or. iblank .eq. 1) go to 52 if(pline(1:1) .eq. ',') then istart = 2 else istart = 1 end if done = .false. 53 continue icomma = index(pline(istart:80),',') + istart - 1 if(icomma.eq. istart-1) then icomma = index(pline(istart:80),' ') + istart - 1 done = .true. end if iprop = iprop + 1 prop(iprop) = pline(istart:icomma-1) write(*,'(a,i3,2x,a,1x,a)')' Property',iprop,'is', . prop(iprop)(1:plen(prop(iprop))) c c Flag which we have do i = 1, MAXPROP if(prop(iprop)(1:plen(prop(iprop))) .eq. . gp(i)(1:plen(gp(i)))) have_prop(i) = .true. end do istart = icomma + 1 if(.not.done) go to 53 52 continue c c c Now assemble a new property list do i = 1, 80 pline(i:i) = ' ' end do c istart = 1 do i = 1, igprop - 1 if(have_prop(i)) then pline(istart:) = ',' istart = istart + 1 pline(istart:) = gp(i)(1:plen(gp(i))) istart = istart + plen(gp(i)) end if end do c pline(istart:) = ',' istart = istart + 1 pline(istart:) = gp(igprop)(1:plen(gp(igprop))) istart = istart + plen(gp(igprop)) c do i = igprop + 1 , MAXPROP if(have_prop(i)) then pline(istart:) = ',' istart = istart + 1 pline(istart:) = gp(i)(1:plen(gp(i))) istart = istart + plen(gp(i)) end if end do c c Write out the updated list of properties c write(*,*) write(*,*) 'Updated property list is:',pline write(igraspout) pline c c Now get the rest c read(igraspin)line read(line,*) nvert,ntri,ngrid write(*,*) 'Number of vertices= ',nvert write(*,*) 'Number of triangles=',ntri write(*,*) 'Size of grid= ',ngrid c Do some bound checking on MAXVERT for nvert and ntri if(nvert.gt.MAXVERT) then write(*,fmt='(a,/,a)') . ' Number of vertices exceeds MAXVERT', . ' Increase parameter MAXVERT and recompile' call ccperr(1, . ' Subroutine merge_with_grasp: nvert exceeds MAXVERT') end if if(3*ntri.gt.MAXVERT) then write(*,fmt='(a,/,a)') . ' Number of triangles exceeds MAXVERT/3', . ' Increase parameter MAXVERT and recompile' call ccperr(1, . ' Subroutine merge_with_grasp: ntri*3 exceeds MAXVERT') end if write(igraspout) line read(igraspin) line write(*,*) 'Midpoint',line(1:plen(line)) write(igraspout) line call read_real_line(nvert*3,vert,igraspin) call write_real_line(nvert*3,vert,igraspout) write(*,*) 'vertices read' call read_real_line(nvert*3,buffer,igraspin) call write_real_line(nvert*3,buffer,igraspout) write(*,*) 'accessibles read' call read_real_line(nvert*3,buffer,igraspin) call write_real_line(nvert*3,buffer,igraspout) write(*,*) 'normals read' if(iform .eq. 2) then call read_i4_line(ntri*3,ibuffer4,igraspin) call write_i4_line(ntri*3,ibuffer4,igraspout) write(*,*) 'triangles read...format=2' else call read_i2_line(ntri*3,ibuffer4,igraspin) call write_i2_line(ntri*3,ibuffer4,igraspout) write(*,*) 'triangles read...format=1' end if c write(*,*) c c c Now in read properties all properties up to gproperty1 c and write them out c do i = 1, igprop - 1 if(have_prop(i)) then write(*,*) 'Reading ',gp(i)(1:plen(gp(i))) call read_real_line(nvert,buffer,igraspin) call write_real_line(nvert,buffer,igraspout) write(*,*) 'Writing ',gp(i)(1:plen(gp(i))) end if end do c c Compute S value and assign to General Property 1 c ------------------------------------------------ c Find nearest point for this surface do l = 1, nvert dmin = 9.e20 do i = 1, m ip = mi(i) if(doflag(ip)) then d = (xp(1,ip) - vert(1,l))**2 + . (xp(2,ip) - vert(2,l))**2 + . (xp(3,ip) - vert(3,l))**2 if(d .lt. dmin) then dmin = d c if you want darray put it in here s = - rn(ip) end if end if end do if(sqrt(dmin) .ge. dmatch) s = smin curv(l) = s end do c c Write out the new general property 1 values if(have_prop(igprop)) then write(*,*) 'Reading ',gp(igprop)(1:plen(gp(igprop))) call read_real_line(nvert,buffer,igraspin) end if write(*,*) 'Writing gproperty with S values' call write_real_line(nvert,curv,igraspout) c c ------------------------------------------------------- c c Are there still values to read from the file? do i = igprop + 1, MAXPROP if(have_prop(i)) then write(*,*) 'Reading ',gp(i)(1:plen(gp(i))) call read_real_line(nvert,buffer,igraspin) write(*,*) 'Writing ',gp(i)(1:plen(gp(i))) call write_real_line(nvert,buffer,igraspout) end if end do read(1,err=100,end=100) c1 write(*,*)' There appear to be additional lines in the GRASP file' write(*,*)' Continuing anyway...' 100 write(*,*) write(*,45)ngrasp 45 format(' Successful merge of grasp file',I2) return end subroutine read_real_line(n,a,iun) real a(n) integer iun read(iun,err=10,end=11) a return 10 call ccperr(1, 'read error in s/r read_real_line') stop 11 write(*,*) 'eof error in s/r read_real_line' stop end subroutine read_i2_line(n,a,iun) integer*2 a(n) integer iun read(iun) a return end subroutine read_i4_line(n,a,iun) integer*4 a(n) integer iun read(iun) a return end subroutine write_real_line(n,a,iun) real a(n) integer iun write(iun) a return end subroutine write_i2_line(n,a,iun) integer*2 a(n) integer iun write(iun) a return end subroutine write_i4_line(n,a,iun) integer*4 a(n) integer iun write(iun) a return end subroutine summary c common/SUMM/Dab,Sab real*8 Dab(4),Sab(4) c write(*,10)Dab(1),Dab(2),(Dab(1)+Dab(2))/2.0,Dab(3),Dab(4), . (Dab(3)+Dab(4))/2.0, . Sab(1),Sab(2),(Sab(1)+Sab(2))/2.0,Sab(3),Sab(4), . (Sab(3)+Sab(4))/2.0,(Sab(3)+Sab(4))/2.0 10 format(//, . ' Summary of results: ',/, . '________________________________________________________',/ . ' D(A->B) D(B->A) D(A->B)+D(B->A)/2 ',/, . ' Mean ',1F8.3,' ',1F8.3,' ', 1F8.3,/, . ' Median ',1F8.3,' ',1F8.3,' ', 1F8.3,//, . ' S(A->B) S(B->A) S(A->B)+S(B->A)/2',/, . ' Mean ',1F8.3,' ',1F8.3,' ', 1F8.3,/, . ' Median ',1F8.3,' ',1F8.3,' ', 1F8.3,/,/, . ' Shape complementarity statistic Sc = ', 1F8.3,/ . '________________________________________________________',/) c return end