program align_frag changes made 18 September 1996 *** c double precision rotation matrix and centers of gravity c introduced a suggestion of T N Bhat to improve reliability c of report of angles when the overall rotation is small c (09-May-1988). c c program to align a subset of each of two coordinate sets using c the Kabsch algorithm. After that, produce a revised list of c matches using the Needleman-Wunsch algorithm. c c option to do both input coordinate sets globally or by a c probe approach. added at Weizmann Institute, July 1984. c c option to utilize sequence information during alignment c added at Weizmann Institute, August 1984. c c add more information in "show" file as requested by M Levitt c write a parallel file having the sequence alignment in a c line by line array, September 1984. c c provide means to "repair" the pair list; perhaps it is c deficient due to freezing in input which differs in length. c will be done when freeze is first command. c c provide means to specify an initial trial alignment from which c the work will proceed instead of using even_up_list_length c c probe and sequence options removed March 1996. c parameter (maxnat = 2500) common rot,cg0,cg1,x0,x1 real*8 rot(3,3),cg0(3),cg1(3) real d(3), 1 x0(3,maxnat),x1(3,maxnat),x2(3,maxnat),x3(3,maxnat),x4(3) real penalty/0.0/,new_penalty/0.0/ integer process_option,auto_cycle,cycle_limit integer*2 mark_pen integer used_1(maxnat),used_2(maxnat) integer trial_strings,to_keep(100,2) character*80 trial_string integer*2 mark_hold integer*2 marks(maxnat,maxnat) integer*2 row_max(maxnat) real scratch(maxnat/2,maxnat/2) equivalence (marks(1,1),scratch(1,1)) character*1 one_character character*1 sequence_1(maxnat),sequence_2(maxnat) character*4 flag character*10 names_1(maxnat),names_2(maxnat) character*20 date_time character*25 plus,stars character*72 line_in character*250 file1,file2,file3,file4,file5 logical any_name,auto,enough,freeze,length_check,loners logical noshow/.true./,nowrite/.true./,pdb,show,trial,turn,write data cycle_limit/25/,smult/3.0/,smax/10.0/,turn/.false./ c initialize rotation matrix to identity do j=1,3 do i=1,3 rot(i,j) = 0.0d0 end do rot(j,j) = 1.0d0 end do call lib$date_time(date_time) write (6,100) '0Align program of G H Cohen, NIH, Bethesda, MD' write (6,100) ' Reference is JMB (1986), 190, 593-604' write (6,100) ' and J Appl Cryst (1997), 30, 1160-1161' write (6,100) ' We greet you on ',date_time 100 format (a,a) 110 format (a,$) process_option = 0 do while (process_option.lt.1 .or. process_option.gt.4) write (6,100) write (6,100) ' You may: 1 - use all atoms' write (6,100) ' 2 - use main-chain atoms' write (6,100) ' 3 - use CA atoms only' write (6,100) ' 4 - use main-chain atoms with CB' write (6,100) ' ' write (6,110) ' Enter option (1, 2 , 3 or 4) ' read (5,240,err=10) process_option write (*,120) ' Option',process_option,' was chosen' 120 format (a,i2,a) 10 continue end do if (process_option.ne.3) then write (6,110) ' Require atom name match ? [Yes] ' read (5,100) one_character call str_upcase(one_character,one_character) any_name = one_character.eq.'N' end if c read coordinate set 1 call read_coordinates(1,file1,process_option,names_1,x2,x4, 1 used_1,n_used_1,i,maxnat) n_read_1 = n_used_1 c read coordinate set 2 call read_coordinates(2,file2,process_option,names_2,x3,x4, 1 used_2,n_used_2,n_hold,maxnat) n_read_2 = n_used_2 write (6,100) ' ' write (6,110) ' Supply trial alignment? [y/N] ' read (5,100) flag trial = flag(1:1).eq.'Y' .or. flag(1:1).eq.'y' .or. 1 flag(1:1).eq.'T' .or. flag(1:1).eq.'t' if (trial) then write (6,100) ' ' write (6,130) 130 format ( 1 ' Specify a range of CA atoms in the fixed molecule followed by'/ 1 ' the first CA atom in the corresponding range for the moving'/ 1' molecule. Specify one range per entry. The range may contain'/ 1 ' a residue number or a chain id followed by a range. If a'/ 1 ' chain id is given, it must be followed immediately by a colon'/ 1 ' and a space [as in - A: 1 100]. If no chain id is specified'/ 1' and the file contains more than one chain the first occurrence'/ 1 ' of a residue having that number will be used.'/ 1 ' A blank line terminates the loop.'/) nc = 1 tk1 = 0 trial_strings = 0 do while (nc.gt.0 .or. nc.eq.-2) write (6,100) ' Enter atom range for the fixed molecule:' read (5,200) nc,trial_string if (nc.gt.0) then if (trial_string.ne.' ') then call read_string(trial_string,i1,i2,names_1,n_read_1) if (i1.gt.0) then tk1 = tk1 + 2 to_keep(tk1-1,1) = i1 if (i2.gt.0 .and. i2.ge.i1) then to_keep(tk1,1) = i2 else to_keep(tk1,1) = i1 end if write (6,100) 1 ' Enter first atom in range for the moving molecule:' read (5,200) nc,trial_string if (nc.gt.0) then if (trial_string.ne.' ') then call read_string(trial_string,i1,i2, 1 names_2,n_read_2) if (i1.gt.0) then to_keep(tk1-1,2) = i1 to_keep(tk1,2) = i1 + 1 (to_keep(tk1,1)-to_keep(tk1-1,1)) else tk1 = tk1 - 2 trial_strings = trial_strings - 1 if (i1.eq.-7 .and. i2.eq.-7) then nc = -2 else nc = -1 end if end if end if end if trial_strings = trial_strings + 1 else if (i1.eq.-7 .and. i2.eq.-7) then nc = -2 else nc = -1 end if end if end if end if end do length_check = .false. n_used_1 = 0 n_used_2 = 0 do i=1,trial_strings do j=to_keep(2*i-1,1),to_keep(2*i,1) if (names_1(j)(6:9).eq.'CA ') then n_used_1 = n_used_1 + 1 used_1(n_used_1) = j end if end do do j=to_keep(2*i-1,2),to_keep(2*i,2) if (names_2(j)(6:9).eq.'CA ') then n_used_2 = n_used_2 + 1 used_2(n_used_2) = j end if end do if (n_used_1.ne.n_used_2) then write (6,140) n_used_1 write (6,140) (used_1(ii),ii=1,n_used_1) write (6,140) n_used_2 write (6,140) (used_2(ii),ii=1,n_used_2) 140 format (x15i5) stop 'List length mismatch' end if end do else if (n_used_1.ne.n_used_2) then length_check = .true. call even_up_list_length(used_1,used_2,n_used_1,n_used_2, 1 maxnat) else length_check = .false. end if end if auto = .false. write = .false. show = .false. loners = .false. do while (.true.) if (.not.auto) then write (6,100) write (6,100) ' The current global values are:' write (6,150) smult,smax 150 format (' Limits: smult,smax=',f6.2,',',f6.2) write (6,160) penalty 160 format (' Penalty =',f5.1) if (loners) then write (6,100) ' Singles are permitted (ON)' else write (6,100) ' Singles are disallowed (OFF)' end if if (show .and. .not.noshow) then write (6,100) ' Name of "SHOW" file is: ', 1 file4(1:index(file4//' ',' ')-1) end if if (write .and. .not.nowrite) then write (6,100) ' Name of rotated coordinate file is: ', 1 file3(1:index(file3//' ',' ')-1) end if flag = ' ' do while (flag.eq.' ') write (6,100) write (6,100) ' Enter option: AUTO, FREEze, GO,' write (6,100) ' HELP, LIMIt, MATRix,' write (6,100) ' PAIRs, PENAlty, REMOve,' write (6,100) ' SHOW, SINGles, STOP,' write (6,110) ' WRITe ? ' read (5,100,end=20) flag end do go to 30 20 flag = 'STOP' 30 call str_upcase(flag,flag) if (length_check .and. flag.ne.'FREE') then length_check = .false. end if else flag = ' ' end if if (flag.eq.'PAIR') then write (6,170) 1 (used_1(i),used_2(i),i=1,max(n_used_1,n_used_2)) 170 format (/' The pairs are:'/(8(i6,i4))) else if (flag.eq.'PENA') then new_penalty = -1.0 write (6,180) penalty 180 format 1 (/' Current penalty is',f6.1,' Ang. Enter new value: '$) do while (new_penalty.lt.0.0) read (5,'(q,f10.0)') nc,new_penalty if (new_penalty.lt.0.0) then write (6,100) ' Penalty must not be negative' write (6,180) penalty end if end do if (nc.gt.0) then if (new_penalty.le.smax) then penalty = new_penalty else write (6,100) ' The penalty may not exceed "smax"' end if end if else if (flag.eq.'HELP') then write (6,190) 190 format (/ 1' To enter an option it is sufficient to type only the ', 1 'capitalized letters.'/ 1 / 1' AUTO - Cycle through alignment procedure until convergence ', 1 '(max 25 cycles)'/ 1' FREEze - Show agreement between two sets of coordinates ', 1 'without moving either'/ 1' GO - Do one cycle of the alignment procedure'/ 1' HELP - Print this message'/ 1' LIMIt - Reset the limits that determine which pairs of ', 1 'potential matches'/ 1' are to be included. Smax is the maximum distance; ', 1 'Smult is the'/ 1' multiplier on the rms agreement beyond which a pair ', 1 'is ignored'/ 1' MATRix - Introduce a rotation matrix for the next refinement ', 1 'cycle'/ 1' PAIRs - Show which pairs are currently included in the ', 1 'alignment'/ 1' PENAlty - Set a gap penalty to make gaps unfavorable'/ 1' REMOve - remove one or more pairs from the current list of ', 1 'pairs'/ 1' SHOW - Set a name for the file which will contain the ', 1 'detailed results of'/ 1' the refinement'/ 1' SINGles - Allow/disallow isolated pairs [strings less than 2 ', 1 'in length]'/ 1' to contribute to the refinement'/ 1' STOP - Exit program'/ 1' WRITe - Set a name for the file to which the rotated moving ', 1 'molecule may be'/ 1' written'/) write (6,110) ' Press return to continue - ' read (5,100) flag else if (flag.eq.'WRIT') then write = .true. write (6,110) 1 ' Enter the name for the rotated coordinate file: ' read (5,200) nc,file3 200 format (q,a) nowrite = nc.eq.0 else if (flag.eq.'SHOW') then show = .true. write (6,110) ' Enter the name for the "SHOW" file: ' read (5,200) nc,file4 noshow = nc.eq.0 else if (flag.eq.'AUTO') then auto_cycle = 1 auto = .true. np_old = -1 sdp_old = 1.e20 write (6,110) 1 ' Enter the name for the rotated coordinate file: ' read (5,200) nc,file3 nowrite = nc.eq.0 write (6,110) ' Enter the name for the "SHOW" file: ' read (5,200) nc,file4 noshow = nc.eq.0 else if (flag.eq.'SING') then if (loners) then flag = 'ON ' write (6,100) ' SINGLES are permitted' else flag = 'OFF ' write (6,100) ' SINGLES are NOT permitted' end if write (6,100) ' SINGLES are now '//flag(1:3) flag = ' ' do while (flag.ne.'ON ' .and. flag.ne.'OFF ') write (6,110) 1 ' Turn SINGLES on or off (respond ON or OFF) ' read (5,100) flag call str_upcase(flag,flag) end do loners = flag.eq.'ON ' else if (flag.eq.'MATR') then do while (.not.turn) write (6,100) ' Enter a rotation matrix (one row/line) -' read (5,210) ((rot(i,j),j=1,3),i=1,3) 210 format (3f10.0) k = 0 do j = 1,3 do i = 1,3 if (rot(i,j).ne.0.0d0) k=1 end do end do turn = k.eq.1 end do else if (flag.eq.'STOP') then call lib$date_time(date_time) write (6,100) ' ' write (6,100) ' We take leave on ',date_time write (6,100) ' ' call exit else if (flag.eq.'LIMI') then write (6,220) ' Replacing current smult,smax=',smult,smax 220 format (a,f6.2,',',f6.2,': ',$) read (5,230) nc,smult1,smax1 230 format (q,2f10.0) if (nc.gt.0) then smult = smult1 smax = smax1 if (penalty.gt.smax) then write (6,100) ' The penalty will be reduced to', 1 ' agree with smax' penalty = smax end if end if else if (flag.eq.'REMO') then write (6,110) ' Remove from pair - to pair: ' read (5,240) i,j 240 format (2i5) if (i.gt.0 .and. j.gt.0 .and. i.le.n_used_1 .and. 1 j.le.n_used_1) then drop = j - i + 1 do k=j+1,n_used_1 used_1(i) = used_1(k) used_2(i) = used_2(k) i = i + 1 end do n_used_1 = n_used_1 - drop n_used_2 = n_used_1 write (6,170) (used_1(i),used_2(i),i=1,n_used_1) else write (6,250) ' Error in input. i,j,n_used_1=', 1 i,j,n_used_1 250 format (a,3i5) end if else if (auto .or. flag.eq.'FREE' .or. flag.eq.'GO ') then if (flag.eq.'FREE') then freeze = .true. write (6,220) ' Replacing current smult,smax=', 1 smult,smax read (5,230) nc,smult1,smax1 if (nc.gt.0) then smult = smult1 smax = smax1 if (penalty.gt.smax) then write (6,100) ' The penalty will be reduced to', 1 ' agree with smax' penalty = smax end if end if write (6,110) ' Enter the name for the "SHOW" file: ' read (5,200) nc,file4 noshow = nc.eq.0 show = .not.noshow else freeze = .false. length_check = .false. end if c now go to it c load the current set of coordinates do i=1,n_used_1 do j=1,3 x0(j,i) = x2(j,used_1(i)) x1(j,i) = x3(j,used_2(i)) end do end do c if turn, a transformation has been loaded via the matrix c command and the LS fitting will not be done in this pass c the center of gravities which are normally calculated in c KABROT will be required below and are calculated here if (turn) then do i=1,3 cg0(i) = 0.0d0 cg1(i) = 0.0d0 end do do i=1,n_used_1 do k=1,3 cg0(k) = cg0(k) + x0(k,i) cg1(k) = cg1(k) + x1(k,i) end do end do do i=1,3 cg0(i) = cg0(i)/float(n_used_1) cg1(i) = cg1(i)/float(n_used_1) end do end if if (.not.(freeze.or.turn)) then c apply kabsch algorithm to determine rotation call kabrot(n_used_1,x0,x1) end if turn = .false. c transform the coordinates which were used c according to the result obtained do i=1,n_used_2 do j=1,3 aaa = cg0(j) do k=1,3 aaa = aaa + rot(j,k)*(x1(k,i)-cg1(k)) end do d(j) = aaa end do do j=1,3 x1(j,i) = d(j) end do end do c and measure the agreement iamx=0 sump=0. np=0 dmx=0. do i=1,n_used_1 do j=1,3 d(j)=x1(j,i)-x0(j,i) end do dd = d(1)*d(1)+d(2)*d(2)+d(3)*d(3) if (dd.gt.dmx) then dmx=dd iamx=i end if sump=sump+dd np=np+1 end do sdp=sump/np sdp=sqrt(sdp) dmx=sqrt(dmx) if (auto) then enough = (sdp.eq.sdp_old .and. np.eq.np_old) if (enough .or. (auto_cycle.ge.cycle_limit)) then show = .true. write= .true. auto = .false. if (auto_cycle.ge.cycle_limit) then write (6,260) cycle_limit 260 format (//' *** WARNING *** WARNING ***'// 1 ' The auto_cycle_limit (',i3, 1 ') has been exhausted.'/ 1 ' If progress is still evident, you ', 1 'may continue. If the problem is'/ 1 ' oscillating, however, you should ', 1 'consider terminating.'/) end if else sdp_old = sdp np_old = np auto_cycle = auto_cycle + 1 end if end if write (6,270) np,sump,sdp,dmx,iamx 270 format (/' Result of alignment: Sum of dev**2 ', 1 'rms dev.' 1 /' For',i5,' target pairs',4x,f10.3,4x,f10.5,:, 1 /' Maximum distance =',f10.5,' at pair',i5) if (enough .or. .not.auto) then call tell_transformation (6) end if c transform all of the moving coordinates do i=1,n_read_2 do j=1,3 aaa = cg0(j) do k=1,3 aaa = aaa + rot(j,k)*(x3(k,i)-cg1(k)) end do x1(j,i) = aaa end do end do if (write .and. .not.nowrite) then open (unit=1,file=file3,status='new', 1 carriagecontrol='list') inquire (unit=1,name=file3) write (6,280) file3(:index(file3//' ',' ')-1) 280 format (/' ',a,' is opened for coordinate output') do i=1,n_hold read (99,100) line_in pdb = line_in(1:6).eq.'ATOM ' if (pdb) then read (line_in(31:54),'(3f8.3)') x4 else read (line_in(16:45),'(3f10.5)') x4 end if do j=1,3 aaa = cg0(j) do k=1,3 aaa = aaa + rot(j,k)*(x4(k)-cg1(k)) end do d(j) = aaa end do if (pdb) then write (line_in(31:54),'(3f8.3)') d else write (line_in(16:45),'(3f10.5)') d end if call str_trim(line_in,line_in//' ',iline_in) write (1,100) line_in(1:iline_in) end do if (pdb) then write (1,100) 'END' end if close (unit=1) rewind 99 write = .false. nowrite = .true. end if c record distances that are within the selection criteria c x1 represents all of x3 after the c last transformation was applied c names_1 goes with x2 (fixed set) c names_2 goes with x3 (or x1) i1 = 1 i2 = n_read_1 j1 = 1 j2 = n_read_2 search_limit = min ( smult*sdp + 0.0001 , smax ) call set_marks(marks,names_1,names_2,x1,x2,i1,i2,j1,j2, 1 search_limit,any_name,maxnat,penalty, 1 mark_pen) if (freeze .and. .not.length_check) then c for pairs so far identified, zero any conflicting c information in the marks array, leaving only the c "frozen" contributions do i=1,n_used_1 i_1 = used_1(i) i_2 = used_2(i) mark_hold = marks(i_1,i_2) do j=i1,i2 marks(j,i_2) = 0 end do do j=j1,j2 marks(i_1,j) = 0 end do marks(i_1,i_2) = mark_hold end do end if if (.not.loners) then call prune_marks(marks,i1,i2,j1,j2,maxnat) end if c first part of Needleman/Wunsch algorithm call sum_marks(marks,row_max,i1,i2,j1,j2,maxnat,mark_pen) c second part of Needleman/Wunsch algorithm, c recording of pairs for the next iteration call trace_marks(marks,i1,i2,j1,j2,used_1,used_2, 1 n_used_1,n_used_2,maxnat,mark_pen) * if (.not.loners .and. (n_used_1.gt.2)) then if (1.eq.2) then ! temporarily disabled c drop lone pairs, even after the pruning c we can be left with lone pairs c we have had problems with the major loop running c off the end of the world, so we trim the ends first do while ((used_1(1).ne.used_1(2)-1 .or. 1 used_2(1).ne.used_2(2)-1) .and. 1 (n_used_1.gt.2)) j = 1 do k=2,n_used_1 used_1(j) = used_1(k) used_2(j) = used_2(k) j = j + 1 end do n_used_1 = n_used_1 - 1 n_used_2 = n_used_1 end do do while ((used_1(n_used_1).ne.used_1(n_used_1-1)+1 .or. 1 used_2(n_used_1).ne.used_2(n_used_1-1)+1) 1 .and. (n_used_1.gt.2)) n_used_1 = n_used_1 - 1 n_used_2 = n_used_1 end do i = 1 do while (i.lt.n_used_1-2 .and. n_used_1.gt.2) i = i + 1 do while ((used_1(i).ne.used_1(i-1)+1 .or. 1 used_2(i).ne.used_2(i-1)+1) .and. 1 (used_1(i).ne.used_1(i+1)-1 .or. 1 used_2(i).ne.used_2(i+1)-1)) j = i do k=i+1,n_used_1 used_1(j) = used_1(k) used_2(j) = used_2(k) j = j + 1 end do n_used_1 = n_used_1 - 1 n_used_2 = n_used_1 end do end do end if if (show .and. .not.noshow) then open (unit=1,file=file4,status='new') inquire (unit=1,name=file4) write (6,290) file4(:index(file4//' ',' ')-1) 290 format (/' ',a,' is opened to show matches') if ((index(file2//' ',' ')+ 1 index(file1//' ',' ')+20).lt.130) then write (1,300) file2(:index(file2//' ',' ')-1), 1 file1(:index(file1//' ',' ')-1) 300 format (' ',a,:,' is fitted to ',:,a) else write (1,300) file2(:index(file2//' ',' ')-1) write (1,300) file1(:index(file1//' ',' ')-1) end if write (1,270) np,sump,sdp call tell_transformation(1) write (1,310) 310 format (/' Fixed Moving Distance'/) do k=1,len(stars)-1 stars(k:k) = '*' end do do k=1,len(plus)-1 plus(k:k) = '+' end do stars(len(stars):len(stars)) = '>' plus(len(plus):len(plus)) = '>' i = min ( len(stars) , int ( sdp/0.2 + 0.5 ) ) stars(i:i) = '|' sump=0. np=0 dmx=0. i1 = 0 j1 = 0 i_seq = 0 l_name = len(names_1(1)) c the following kludge is used to handle the tail c of the molecule in the same framework as the rest used_1(n_used_1+1) = n_read_1 + 1 used_2(n_used_2+1) = n_read_2 + 1 do k=1,n_used_1+1 i = used_1(k) j = used_2(k) if (i.ne.i1+1 .or. j.ne.j1+1) then i0 = i - i1 - 1 j0 = j - j1 - 1 if (i0.gt.0 .and. j0.gt.0) then c calculate all distances between c pairs in the string call get_scratch_distances(scratch,maxnat/2, 1 x2(1,i1+1),names_1(i1+1),i0, 1 x1(1,j1+1),names_2(j1+1),j0, 1 any_name) c for now, the ultimate string length is c limited by the length of the shorter string k0 = min(i0,j0) if (i0.eq.j0) then c for now, if both strings are c equal in length, use as is do k2=1,k0 scratch(k2,1) = scratch(k2,k2) end do call list_pair(names_1(i1+1),names_2(j1+1), 1 sequence_1(i_seq+1),sequence_2(i_seq+1), 1 scratch,'-',k0,plus,process_option) i_seq = i_seq + k0 else c we have to find out which part of x1 or x2 c matches the other string best, c taking k0 at a time small_sum = 9.0e20 is = 0 js = 0 do j2=1,j0-k0+1 do i2=1,i0-k0+1 sum = 0.0 do k2=0,k0-1 sum = sum + scratch(i2+k2,j2+k2) end do if (sum.lt.small_sum) then small_sum = sum is = i2 js = j2 end if end do end do c is the result decent? if (is.eq.0 .or. js.eq.0) 1 stop 'program error with is,js' c now do something about it c report tails that precede the string if (is.gt.1) then i2 = is-1 call list_left_member(names_1(i1+1), 1 sequence_1(i_seq+1), 1 sequence_2(i_seq+1), 1 i2,process_option) i_seq = i_seq + i2 end if if (js.gt.1) then j2 = js-1 call list_right_member(names_2(j1+1), 1 sequence_1(i_seq+1), 1 sequence_2(i_seq+1), 1 j2,l_name,process_option) i_seq = i_seq + j2 end if c now report the forced match do k2=1,k0 scratch(k2,1) = scratch(is+k2-1,js+k2-1) end do call list_pair(names_1(i1+is),names_2(j1+js), 1 sequence_1(i_seq+1),sequence_2(i_seq+1), 1 scratch,'-',k0,plus,process_option) i_seq = i_seq + k0 c now report trailing tails if (i0.gt.is+k0-1) then i2 = i0-is-k0+1 call list_left_member(names_1(i1+is+k0), 1 sequence_1(i_seq+1), 1 sequence_2(i_seq+1), 1 i2,process_option) i_seq = i_seq + i2 end if if (j0.gt.js+k0-1) then j2 = j0-js-k0+1 call list_right_member(names_2(j1+js+k0), 1 sequence_1(i_seq+1), 1 sequence_2(i_seq+1), 1 j2,l_name, 1 process_option) i_seq = i_seq + j2 end if end if else c just show an unmatched segment if (i0.gt.0) then call list_left_member(names_1(i1+1), 1 sequence_1(i_seq+1), 1 sequence_2(i_seq+1), 1 i0,process_option) i_seq = i_seq + i0 end if if (j0.gt.0) then call list_right_member(names_2(j1+1), 1 sequence_1(i_seq+1), 1 sequence_2(i_seq+1), 1 j0,l_name,process_option) i_seq = i_seq + j0 end if end if end if if (k.le.n_used_1) then dx = x2(1,i)-x1(1,j) dy = x2(2,i)-x1(2,j) dz = x2(3,i)-x1(3,j) dd = dx**2 + dy**2 + dz**2 if (dd.gt.dmx) dmx = dd sump = sump + dd i_seq = i_seq + 1 call list_pair(names_1(i),names_2(j), 1 sequence_1(i_seq),sequence_2(i_seq), 1 dd,'=',1,stars,process_option) i1 = i j1 = j end if end do dmx = sqrt(dmx) sdp = sqrt(sump/n_used_1) write (1,320) n_used_1,dmx,sump,sdp 320 format (/' For the',i5,' pairs listed, the maximum ', 1 'distance =',f10.5/' Sum of dev**2 ='f10.3, 1 ', rms deviation ='f10.5) close (unit=1) if (process_option.eq.3) then ! restrict for Ca matching only open (unit=1,file='.seq',status='new', 1 defaultfile=file4,carriagecontrol='list') inquire (unit=1,name=file5) write (6,330) file5(:index(file5//' ',' ')-1) 330 format (' ',a,' is opened to show sequences') write (1,100) file1(:index(file1//' ',' ')-1) 340 format ('Sequence from file: ',a) write (1,350) (sequence_1(i),i=1,i_seq) write (1,100) file2(:index(file2//' ',' ')-1) write (1,350) (sequence_2(j),j=1,i_seq) 350 format (30a2) close (unit=1) end if show = .false. noshow = .true. length_check = .false. end if end if end do end subroutine read_string(string,i1,i2,names,n_read) changes made 18 September 1996 *** implicit none character*(*) string character*(*) names(*) integer n_read character*4 res1,res2 character*1 chain integer lres1,lres2 integer i1,i2,len_str,hyph_pos,i,j,k logical is_chain call str_trim(string,string//' ',len_str) call str_upcase(string,string(1:len_str)) do i=1,len_str if (.not.(ichar(string(i:i)).ge.ichar('0').and. 1 ichar(string(i:i)).le.ichar('9').or. 1 ichar(string(i:i)).ge.ichar('A').and. 1 ichar(string(i:i)).le.ichar('Z').or. 1 string(i:i).eq.':')) then string(i:i) = ' ' end if end do c a chain name will consist of only one non-blank character and c be followed immediately by a colon and at least one blank hyph_pos = index(string(1:len_str)//': ',': ') is_chain = hyph_pos.lt.len_str is_chain = is_chain.and.string(hyph_pos-1:hyph_pos-1).ne.' ' is_chain = is_chain.and. 1 (hyph_pos.eq.2.or.string(hyph_pos-2:hyph_pos-2).eq.' ') if (is_chain) then chain = string(hyph_pos-1:hyph_pos-1) i = hyph_pos + 2 else chain = ' ' i = 1 end if do j = i,len_str if (string(j:j).eq.':') string(j:j) = ' ' end do do while (i.le.len_str .and. string(i:i).eq.' ') i = i + 1 end do if (i.gt.len_str) then i1 = 0 i2 = 0 return end if res1 = ' ' lres1 = 0 do while (i.le.len_str .and. string(i:i).ne.' ') lres1 = lres1 + 1 res1(lres1:lres1) = string(i:i) i = i + 1 end do k = 4 do j = lres1,1,-1 res1(k:k) = res1(j:j) k = k - 1 end do do j = k,1,-1 res1(j:j) = ' ' end do do while (i.le.len_str .and. string(i:i).eq.' ') i = i + 1 end do if (i.gt.len_str) then lres2 = 0 res2 = ' ' end if res2 = ' ' lres2 = 0 do while (i.le.len_str .and. string(i:i).ne.' ') lres2 = lres2 + 1 res2(lres2:lres2) = string(i:i) i = i + 1 end do if (lres2.gt.0) then k = 4 do j = lres2,1,-1 res2(k:k) = res2(j:j) k = k - 1 end do do j = k,1,-1 res2(j:j) = ' ' end do end if i = 1 if (is_chain) then do while (i.lt.n_read .and. chain.ne.names(i)(10:10)) i = i + 1 end do end if do while (i.lt.n_read .and. ((res1.ne.names(i)(2:5) .and. 1 res1(2:4)//' '.ne.names(i)(2:5)) .or. names(i)(6:9).ne.'CA ')) i = i + 1 end do if (names(i)(6:9).eq.'CA ' .and. (res1.eq.names(i)(2:5) .or. 1 res1(2:4)//' '.eq.names(i)(2:5))) then if (is_chain) then if (chain.eq.names(i)(10:10)) then i1 = i else i1 = 0 end if else i1 = i end if else i1 = 0 end if if (i1.ne.0) then write (6,100) 1 'Range starts with '//names(i1)(10:10)//' '//names(i1)(1:9) 100 format (xa) else write (6,100) 'Range start not found in atom list' i1 = -7 i2 = -7 return end if if (lres2.gt.0) then do while (i.lt.n_read .and. ((res2.ne.names(i)(2:5) .and. 1 res2(2:4)//' '.ne.names(i)(2:5)) .or. names(i)(6:9).ne.'CA ')) i = i + 1 end do if (names(i)(6:9).eq.'CA ' .and. (res2.eq.names(i)(2:5) .or. 1 res2(2:4)//' '.eq.names(i)(2:5))) then if (is_chain) then if (chain.eq.names(i)(10:10)) then i2 = i else i2 = 0 end if else i2 = i end if else i2 = 0 end if if (i2.ne.0) then write (6,100) 1 'Range ends with '//names(i2)(10:10)//' '//names(i2)(1:9) else write (6,100) 'Range end not found in atom list' i1 = -7 i2 = -7 return end if end if return end subroutine even_up_list_length(used_1,used_2,n_used_1,n_used_2, 1 maxnat) changes made 29 March 1996 *** integer used_1(maxnat),used_2(maxnat) c be sure that n_used_1 is the same as n_used_2. c if not, drop candidates from the longer list. do while (n_used_1.ne.n_used_2) i = min(n_used_1,n_used_2) j = max(n_used_1,n_used_2) kd = j/(j-i) if (kd.lt.2) kd = 2 k = -kd/2 l = 0 do while (k.le.j .and. l.lt.j-i) k = k + kd l = l + 1 if (n_used_1.gt.n_used_2) then used_1(k) = 0 else used_2(k) = 0 end if end do i = 0 if (n_used_1.gt.n_used_2) then do j=1,n_used_1 if (used_1(j).gt.0) then i = i + 1 used_1(i) = used_1(j) end if end do n_used_1 = i else do j=1,n_used_2 if (used_2(j).gt.0) then i = i + 1 used_2(i) = used_2(j) end if end do n_used_2 = i end if end do return end subroutine read_coordinates(number,file,process_option,names,x,x1, 1 used,n_used,n_hold,maxnat) changes made 29 March 1996 *** integer number,n_used,process_option,used(maxnat) integer this,which_l(2),file_l real d(3),x(3,maxnat),x1(3) character*3 chr_3,get_name character*1 get_chr_1 character*2 chno character*4 main_atoms(5) character*4 atom character*5 residue character*6 which(2) character*72 line_in character*(*) names(maxnat) character*(*) file logical use,pdb data main_atoms/'N','CA','C','O','CB'/ data which/'fixed','moving'/,which_l/5,6/ 100 format (a) 110 format (a,$) 10 write (6,120) which(number)(1:which_l(number)),number 120 format (/' Enter the name of the ',a,' coordinate set (#',i1,') ') write (6,100) 1 ' [The file is assumed to be in PDB format unless it has ]' write (6,100) 1 ' [a : appended to the file name to indicate WAH format. ]' write (6,110) ' File name - ' read (5,100) file call str_trim(file,file,file_l) pdb = .not.file(file_l:file_l).eq.':' if (pdb) then write (6,100) ' Input file will be read as a PDB file' else file(file_l:file_l) = ' ' write (6,100) ' Input file will be read as a WAH file' end if open (unit=1,status='old',readonly,file=file,err=10) write (6,100) ' Enter cell parameters (a,b,c,ca,cb,cg ) or,' write (6,110) ' if coordinates are already cartesian: ' read (5,130) nc,acell,bcell,ccell,ca,cb,cg 130 format (q,6f10.0) if (nc.ne.0) then if (abs(ca).gt.1.0) then write (6,140) ca 140 format (' Input angle of',f8.3,' degrees is converted to' 1 ' the cosine') ca = cosd(ca) end if if (abs(cb).gt.1.0) then write (6,140) cb cb = cosd(cb) end if if (abs(cg).gt.1.0) then write (6,140) cg cg = cosd(cg) end if write (6,150) acell,bcell,ccell,ca,cb,cg 150 format (' cell parameters (a ,b ,c ) =',3f10.5/ 1 ' (ca,cb,cg) =',3f10.5) write (6,*) 'Effective matrix:' sg = sqrt(1.0-cg*cg) cartzz = sqrt(1.0-ca*ca-cb*cb-cg*cg+2.0*ca*cb*cg)/sg cartyz = (cg-cb*ca)/sg geom11 = acell geom12 = bcell*cg geom13 = ccell*cb geom22 = bcell*sg geom23 = ccell*cartyz geom33 = ccell*cartzz write (6,*) geom11, geom12, geom13 write (6,*) 0.0, geom22, geom23 write (6,*) 0.0, 0.0, geom33 end if inquire (unit=1,name=file) write (6,160) file(:index(file//' ',' ')-1),number 160 format (/' ',a,' is opened for coordinate set',i2,' input') if (number.ne.1) then open (unit=99,status='scratch') end if i = 0 ii = 1 n_used = 0 do while (i.le.maxnat) c ".le." permits termination on eof if found in time if (pdb) then line_in = ' ' do while (line_in(1:6).ne.'ATOM ') read (1,100,end=20) line_in end do residue = ' ' read (line_in,170) atom,chr_3,chno,residue(2:5),d 170 format (13x,a,a,a,x,a,3x,3f8.3) else read (1,100,end=20) line_in read (line_in,180) chno,chr_3,residue(1:4),atom,d 180 format (a,x,a,x,2a,3f10.5) end if if (residue(1:1).eq.' ') then if (chr_3.ne.' ') then residue(1:1) = get_chr_1(chr_3) end if end if if (process_option.eq.1) then use = .true. else if (process_option.eq.2) then use = .false. do j=1,4 use = use .or. main_atoms(j).eq.atom end do else if (process_option.eq.3) then use = atom.eq.'CA ' else if (process_option.eq.4) then use = .false. do j=1,5 use = use .or. main_atoms(j).eq.atom end do end if if (use .or. number.ne.1) then if (nc.gt.0) then x1(1) = d(1)*geom11 + d(2)*geom12 + d(3)*geom13 x1(2) = d(2)*geom22 + d(3)*geom23 x1(3) = d(3)*geom33 else x1(1) = d(1) x1(2) = d(2) x1(3) = d(3) end if if (use) then i = i + 1 x(1,i) = x1(1) x(2,i) = x1(2) x(3,i) = x1(3) names(i) = residue//atom//chno(2:2) n_used = n_used + 1 used(n_used) = i end if if (number.ne.1) then if (ichar(residue(2:2)).gt.ichar('9')) then residue(1:1) = residue(2:2) residue(2:2) = ' ' end if if (ichar(residue(3:3)).gt.ichar('9')) then residue(1:1) = residue(3:3) residue(3:3) = ' ' end if if (.not.pdb) then if (chno.eq.' ') then this = ichar(residue(1:1)) if (this.gt.ichar('0') .and. this.le.ichar('9')) 1 then chno = ' '//char(this - ichar('0') + 1) else chno = ' 1' end if end if end if if (chr_3.eq.' ') then chr_3 = get_name(residue(1:1)) end if if (pdb) then line_in(18:20) = chr_3 else line_in(4:6) = chr_3 line_in(8:11) = residue(1:4) end if write (99,100) line_in ii = ii + 1 end if end if end do write (6,190) maxnat 190 format (/' The number of atoms read exceeds the available', 1 ' storage which is',i7/) stop 'Too many atoms' 20 close (unit=1) if (number.eq.1) then write (6,200) number,n_used 200 format (/' For coordinate set',i2,', use',i5,:,', hold',i5) else n_hold = ii - 1 write (6,200) number,n_used,n_hold rewind 99 end if return end character*3 function get_name(short_name) character*(*) short_name character*3 resids(26) data resids/'ALA',' ','CYS','ASP','GLU','PHE','GLY','HIS','ILE', 1 'PCA','LYS','LEU','MET','ASN','WAT','PRO','GLN','ARG', 2 'SER','THR','SUL','VAL','TRP','HEM','TYR',' '/ index = ichar(short_name(1:1)) - ichar('A') + 1 if (index.le.0 .or. index.gt.26) then get_name = ' ' else get_name = resids(index) end if return end character*1 function get_chr_1(chr_3) character*(*) chr_3 character*3 resids(26) data resids/'ALA',' ','CYS','ASP','GLU','PHE','GLY','HIS','ILE', 1 'PCA','LYS','LEU','MET','ASN','WAT','PRO','GLN','ARG', 2 'SER','THR','SUL','VAL','TRP','HEM','TYR',' '/ i = 1 do while (i.le.26) if (chr_3.eq.resids(i)) then get_chr_1 = char(i+ichar('A')-1) i = 27 else i = i + 1 end if end do return end subroutine set_marks(marks,names_1,names_2,x1,x2,i1,i2,j1,j2, 1 search_limit,any_name,maxnat,penalty, 1 mark_pen) changes made 29 March 1996 *** real x1(3,maxnat),x2(3,maxnat) integer*2 marks(maxnat,maxnat) real mark_scale,penalty integer*2 mark_pen character*4 test_atom_2 character*(*) names_1(maxnat),names_2(maxnat) logical any_name i_span = i2 - i1 + 1 j_span = j2 - j1 + 1 search_limit_squared = search_limit**2 mark_scale = 32767./min(i_span,j_span)/search_limit_squared do j=j1,j2 test_atom_2 = names_2(j)(6:9) do i=i1,i2 marks(i,j) = 0 if (any_name.or.(test_atom_2.eq.names_1(i)(6:9))) then dx = abs(x2(1,i)-x1(1,j)) if (dx.le.search_limit) then dy = abs(x2(2,i)-x1(2,j)) if (dy.le.search_limit) then dz = abs(x2(3,i)-x1(3,j)) if (dz.le.search_limit) then dd = dx**2 + dy**2 + dz**2 if (dd.le.search_limit_squared) then dd = sqrt(dd) marks(i,j) = mark_scale * 1 (search_limit-dd)**2 end if end if end if end if end if end do end do mark_pen = mark_scale*penalty**2 return end subroutine prune_marks(marks,i1,i2,j1,j2,maxnat) changes made 23-Jan-1985 *** integer*2 marks(maxnat,maxnat) c prune the marks array to remove any potential marks c which are not part of a chain do j=j1+1,j2-1 do i=i1+1,i2-1 if (marks(i-1,j-1).eq.0 .and. 1 marks(i+1,j+1).eq.0) marks(i,j) = 0 end do end do do i=i1,i2-1 if (marks(i+1,j1+1).eq.0) marks(i,j1) = 0 end do do i=i1+1,i2 if (marks(i-1,j2-1).eq.0) marks(i,j2) = 0 end do do j=j1,j2-1 if (marks(i1+1,j+1).eq.0) marks(i1,j) = 0 end do do j=j1+1,j2 if (marks(i2-1,j-1).eq.0) marks(i2,j) = 0 end do marks(i1,j2) = 0 marks(i2,j1) = 0 return end subroutine sum_marks(marks,row_max,i1,i2,j1,j2,maxnat,mark_pen) changes made 29 March 1996 *** integer*2 col_max,mark_test,max_here integer*2 marks(maxnat,maxnat) integer*2 row_max(maxnat) integer*2 mark_pen do i=i1,i2 row_max(i) = 0 end do do j=j2,j1+1,-1 col_max = 0 do i=i2,i1+1,-1 mark_test = marks(i,j) if (mark_test.gt.row_max(i)) then row_max(i) = mark_test end if if (mark_test.gt.col_max) then col_max = mark_test end if max_here = max(row_max(i),col_max) marks(i-1,j-1) = marks(i-1,j-1) + 1 max(marks(i,j),max_here-mark_pen) end do end do return end subroutine trace_marks(marks,i1,i2,j1,j2,used_1,used_2, 1 n_used_1,n_used_2,maxnat,mark_pen) changes made 2 April 1996 *** integer used_1(maxnat),used_2(maxnat) integer*2 last_maximum,mark_test,maximum_found,mark_pen integer*2 marks(maxnat,maxnat) integer d,shortest n_used_1 = 0 i = i1 j = j1 last_maximum = -1 do while (i.le.i2 .and. j.le.j2) i_best = i j_best = j n_found = 1 maximum_found = marks(i,j) shortest = 0 do k=i+1,i2 mark_test = marks(k,j) - mark_pen if (mark_test.ge.maximum_found) then d = k - i if (d.lt.shortest .or. 1 mark_test.gt.maximum_found .or. 1 mark_pen.eq.0) then if (mark_test.eq.maximum_found .and. 1 mark_pen.eq.0) then n_found = n_found + 1 else n_found = 1 maximum_found = mark_test i_best = k j_best = j shortest = d end if end if end if end do do l=j+1,j2 mark_test = marks(i,l) - mark_pen if (mark_test.ge.maximum_found) then d = l - j if (d.lt.shortest .or. 1 mark_test.gt.maximum_found .or. 1 mark_pen.eq.0) then if (mark_test.eq.maximum_found .and. 1 mark_pen.eq.0) then n_found = n_found + 1 else n_found = 1 maximum_found = mark_test i_best = i j_best = l shortest = d end if end if end if end do if (maximum_found.eq.0) then c force exit from the "do while" i = i2 + 1 else c assume that a multiple find implies a null (this c assumption requires sufficient sensitivity in the c weights) if (n_found.gt.1) then i = i + 1 j = j + 1 else c a repeat maximum implies no match c in the previous position if (n_used_1.eq.0 .or. 1 (maximum_found.ne.last_maximum)) then c this one is a real point to record n_used_1 = n_used_1 + 1 end if c record the currently found position c (or update the current position pointer) used_1(n_used_1) = i_best used_2(n_used_1) = j_best i = i_best + 1 j = j_best + 1 last_maximum = maximum_found end if end if end do n_used_2 = n_used_1 return end subroutine tell_transformation(lun) changes made 3-Aug-2000 *** implicit real*8 (a-h,o-z) common rot,cg0,cg1,separation real*8 dt(3),o(3),r(3),rot(3,3),cg0(3),cg1(3),sp(3),zero/1.0d-6/ save dt,o,r,sp save deltadeg,phideg,thetadeg,psideg,s degtorad = 180.0d0/acos(-1.0d0) if (lun.eq.6) then c what is the original separation of the molecules separation = 0.0 do j=1,3 separation = separation + (cg1(j)-cg0(j))**2 end do separation = sqrt(separation) c consolidate the translational part do j=1,3 dt(j) = cg0(j) do k=1,3 dt(j) = dt(j) - rot(j,k)*cg1(k) end do end do c rephrase the transformation if (abs(rot(3,1)).gt.zero .or. abs(rot(3,2)).gt.zero) then phi = atan2(rot(3,1),-rot(3,2)) else phi = 0.0d0 end if if (abs(rot(1,3)).gt.zero .or. abs(rot(2,3)).gt.zero) then psi = atan2(rot(1,3),rot(2,3)) else psi = 0.0d0 end if if (abs(rot(3,2)).ge.abs(rot(2,3))) then sthet = - rot(3,2)/cos(phi) else sthet = rot(2,3)/cos(psi) end if theta = atan2(sthet,rot(3,3)) sin_h_theta = sin(0.5d0*theta) cos_h_theta = cos(0.5d0*theta) if (theta.ne.0.0d0) then cot_h_theta = cos_h_theta/sin_h_theta else cot_h_theta = 0.0d0 end if sin_h_phipsi = sin(0.5d0*(phi+psi)) sin_h_phimpsi = sin(0.5d0*(phi-psi)) cos_h_phimpsi = cos(0.5d0*(phi-psi)) q = sqrt(1.0d0+(cot_h_theta*sin_h_phipsi)**2) if (q.ne.0.0d0) then r(1) = cos_h_phimpsi/q r(2) = sin_h_phimpsi/q r(3) = cot_h_theta*sin_h_phipsi/q else r(1) = 0.0d0 r(2) = 0.0d0 r(3) = 0.0d0 end if arg = 0.5d0*(rot(1,1)+rot(2,2)+rot(3,3)-1.0d0) arg = sign(min(abs(arg),1.0d0),arg) cos_h_delta = cos(0.5d0*acos(arg)) sin_h_delta = r(1) * sin_h_theta * cos_h_phimpsi 1 + r(2) * sin_h_theta * sin_h_phimpsi 1 + r(3) * cos_h_theta * sin_h_phipsi delta = 2.0d0 * atan2(sin_h_delta,cos_h_delta) phideg = phi*degtorad thetadeg = theta*degtorad psideg = psi*degtorad deltadeg = delta*degtorad if (delta.ne.0.0d0) then cot_h_delta = cos_h_delta/sin_h_delta else cot_h_delta = 0.0d0 end if c translation along axis s = r(1)*dt(1) + r(2)*dt(2) + r(3)*dt(3) c translation perpendicular to axis do i=1,3 sp(i) = cg1(i) - cg0(i) - s * r(i) end do c point on axis (McLachlan, JMB(1979) 128, 49-79) do i=1,3 j = mod(i,3) + 1 k = mod(i+1,3) + 1 o(i) = 0.5d0 * ( ( r(j) * sp(k) - r(k) * sp(j) ) 1 * cot_h_delta + cg1(i) + cg0(i) ) end do end if write (lun,99) separation 99 format (/' The centroids of the input coordinate sets', 1 ' are separated by:',f9.5,' Angstroms') write (lun,100) cg0,cg1,((rot(i,j),j=1,3),i=1,3) 100 format (/' Coordinates of centroids:'/ 1 /' Fixed'f10.5,2f11.5,' Moving'f10.5,2f11.5/ 1 /' Rotation matrix applied to latest set:' 1 6x3f11.6,2(//45x3f11.6)) write (lun,110) dt 110 format (/' The combined translation vector is:'3f11.5) write (lun,120) deltadeg,r,phideg,thetadeg,psideg 120 format (/' The matrix corresponds to a rotation of'f10.5, 1 ' degrees' 1 /' around the direction'3f10.5/ 1 /' The equivalent Eulerean rotation is' 1 /' Phi =',f10.5,', Theta =',f10.5,', Psi ='f10.5/) write (lun,130) o,s 130 format (' A point on the rotation axis is'3f10.5/ 1 ' The moving model translates'f10.5, 1 ' along the axis during rotation') return end subroutine list_pair(names_1,names_2,sequence_1,sequence_2,dis, 1 sign,n,stars,process_option) changes made 29 March 1996 *** character*(*) names_1(n),names_2(n),sequence_1(n),sequence_2(n), 1 sign,stars character*1 one_character real dis(n) integer process_option do i=1,n dd = sqrt(dis(i)) n_stars = min(int(dd/0.2+0.5),len(stars)) if (names_1(i)(1:1).eq.names_2(i)(1:1)) then one_character = sign else one_character = ' ' end if write (1,100) names_1(i)(10:10),names_1(i)(1:9),one_character, 1 names_2(i)(10:10),names_2(i)(1:9), 1 dd,stars(1:n_stars) 100 format (1x,a,x,a,a,1x,a,x,a,f7.2,1x,a) if (process_option.eq.3) then if (names_1(i)(6:7).eq.'CA' .and. 1 names_2(i)(6:7).eq.'CA') then sequence_1(i) = names_1(i)(1:1) sequence_2(i) = names_2(i)(1:1) end if end if end do end subroutine list_left_member(names_1,sequence_1,sequence_2,n, 1 process_option) changes made 29 March 1996 *** character*(*) names_1(n),sequence_1(n),sequence_2(n) character*8 name_out integer process_option do i=1,n call str_trim(name_out,names_1(i)(1:9),k) write (1,100) names_1(i)(10:10),name_out(1:k) 100 format (1x,a,x,a) if (process_option .eq. 3) then if (names_1(i)(6:7).eq.'CA') then sequence_1(i) = names_1(i)(1:1) sequence_2(i) = '-' end if end if end do return end subroutine list_right_member(names_2,sequence_1,sequence_2,n, 1 l_name,process_option) changes made 29 March 1996 *** character*(*) names_2(n),sequence_1(n),sequence_2(n) character*8 name_out integer process_option do j=1,n call str_trim(name_out,names_2(j)(1:9),k) write (1,100) names_2(j)(10:10),name_out(1:k) 100 format (x,a,x,a) if (process_option.eq.3) then if (names_2(j)(6:7).eq.'CA') then sequence_1(j) = '-' sequence_2(j) = names_2(j)(1:1) end if end if end do return end subroutine get_scratch_distances(scratch,n_s,x,xn,n_x,y,yn,n_y, 1 any_name) changes made 29 March 1996 *** character*(*) xn(n_x),yn(n_y) character*4 test_name real scratch(n_s,n_s),x(3,n_x),y(3,n_y) logical any_name do j=1,n_y test_name = yn(j)(6:9) do i=1,n_x if (any_name .or. (test_name.eq.xn(i)(6:9))) then scratch(i,j) = (x(1,i)-y(1,j))**2 + 1 (x(2,i)-y(2,j))**2 + 1 (x(3,i)-y(3,j))**2 else scratch(i,j) = 999.0 end if end do end do return end subroutine kabrot(nat,x0,x1) changes made 09-May-1988 *** c c calculates orthogonal rotation matrix using the algorithm of c Kabsch ACTA A32, 922 (1976) c ACTA A34, 827 (1978) c implicit integer *4 (i-n), real*8 (a-h,o-z) common rot,cg0,cg1 real*8 rot(3,3),cg0(3),cg1(3) real*4 x0(3,nat),x1(3,nat) real*8 u(3,3),r(3,3),rtr(3,3),eig(3,3),eigval(3),rtrcol(6),b(3,3) real*8 a3(3),b3(3) c calculate centroids do i=1,3 cg0(i) = 0.0d0 cg1(i) = 0.0d0 end do do i=1,nat do k=1,3 cg0(k) = cg0(k) + x0(k,i) cg1(k) = cg1(k) + x1(k,i) end do end do do i=1,3 cg0(i) = cg0(i)/float(nat) cg1(i) = cg1(i)/float(nat) end do c set up r rotation matrix do i=1,3 do j=1,3 r(i,j) = 0.0d0 do k=1,nat r(i,j) = (x0(i,k)-cg0(i)) * 1 (x1(j,k)-cg1(j)) + r(i,j) end do end do end do c now for rtr (r-transpose * r) do i=1,3 do j=1,3 rtr(i,j)=0.0d0 do k=1,3 rtr(i,j)=rtr(i,j)+r(k,i)*r(k,j) end do end do end do c find eigenvalues and eigenvectors k = 0 do i=1,3 do j=1,i k = k + 1 rtrcol(k) = rtr(i,j) end do end do call eigen(rtrcol,eig,3,0) eigval(1) = rtrcol(1) eigval(2) = rtrcol(3) eigval(3) = rtrcol(6) c eigenvalues and vectors have been sorted by decreasing c size in subroutine EIGEN enabling treatment of degenerate c and zero eigenvalues call cross(eig) if(eigval(3).le.0.0d0) eigval(3)=1.0d0 c calculate Kabsch`s B matrix do i=1,3 do j=1,3 b(j,i)=0.0d0 do k=1,3 b(j,i)=b(j,i)+r(j,k)*eig(k,i) end do b(j,i)=b(j,i)/sqrt(eigval(i)) end do end do c treat degenerate and zero eigenvalues (disabled because c the cross seems to handle the complete problem) c b3(1)=b(1,3) c b3(2)=b(2,3) c b3(3)=b(3,3) call norm(b) call cross(b) c a3(1)=b(1,3) c a3(2)=b(2,3) c a3(3)=b(3,3) c arg=a3(1)*b3(1)+a3(2)*b3(2)+a3(3)*b3(3) c if(arg.le.0.0) then c b(1,3)=-b(1,3) c b(2,3)=-b(2,3) c b(3,3)=-b(3,3) c end if c calculate the u matrix do i=1,3 do j=1,3 u(i,j)=0.0d0 do k=1,3 u(i,j)=b(i,k)*eig(j,k)+u(i,j) end do end do end do det = u(1,1)*u(2,2)*u(3,3) + u(2,1)*u(3,2)*u(1,3) + 1 u(3,1)*u(1,2)*u(2,3) - u(3,1)*u(2,2)*u(1,3) - 1 u(2,1)*u(1,2)*u(3,3) - u(1,1)*u(3,2)*u(2,3) if (det.lt.0.0d0) then write (6,*) ' *** WARNING *** - The matrix is left handed' write (6,*) ' *** WARNING *** - Determinant =',det end if c copy the matrix into rot to return it do j=1,3 do i=1,3 rot(i,j) = u(i,j) end do end do return end subroutine eigen(a,r,n,mv) changes made 09-May-1988 *** implicit real*8 (a-h,o-z), integer*4 (i-n) real*8 a(1),r(1) RANGE=1.0D-6 IF(MV-1) 10,40,10 10 IQ=-N DO 30 J=1,N IQ=IQ+N DO 30 I=1,N IJ=IQ+I R(IJ)=0.0D0 IF(I-J) 30,20,30 20 R(IJ)=1.0D0 30 CONTINUE 40 ANORM=0.0D0 DO 60 I=1,N DO 60 J=I,N IF(I-J) 50,60,50 50 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 60 CONTINUE IF(ANORM) 320,320,70 70 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) IND=0 THR=ANORM 80 THR=THR/FLOAT(N) 90 L=1 100 M=L+1 110 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ IF( ABS(A(LM))-THR) 250,120,120 120 IND=1 LL=L+LQ MM=M+MQ X=0.5D0*(A(LL)-A(MM)) Y=-A(LM)/ SQRT(A(LM)*A(LM)+X*X) IF(X) 130,140,140 130 Y=-Y 140 SINX=Y/ SQRT(2.0D0*(1.0D0+( SQRT(1.0D0-Y*Y)))) SINX2=SINX*SINX COSX= SQRT(1.0D0-SINX2) COSX2=COSX*COSX SINCS =SINX*COSX ILQ=N*(L-1) IMQ=N*(M-1) DO 240 I=1,N IQ=(I*I-I)/2 IF(I-L) 150,220,150 150 IF(I-M) 160,220,170 160 IM=I+MQ GO TO 180 170 IM=M+IQ 180 IF(I-L) 190,200,200 190 IL=I+LQ GO TO 210 200 IL=L+IQ 210 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 220 IF(MV-1) 230,240,230 230 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 240 CONTINUE X=2.0D0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 250 IF(M-N) 260,270,260 260 M=M+1 GO TO 110 270 IF(L-(N-1)) 280,290,280 280 L=L+1 GO TO 100 290 IF(IND-1) 310,300,310 300 IND=0 GO TO 90 c c compare threshold with final norm c 310 IF(THR-ANRMX) 320,320,80 c c sort eigenvalues and eigenvectors c 320 IQ=-N DO 360 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 360 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 330,360,360 330 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 340,360,340 340 DO 350 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 350 R(IMR)=X 360 CONTINUE RETURN END SUBROUTINE CROSS(A) changes made 09-May-1988 *** C C 3RD COLUMN OF A IS CROSS PRODUCT OF FIRST TWO COLUMNS C REAL*8 A(3,3) A(1,3)=A(2,1)*A(3,2)-A(2,2)*A(3,1) A(2,3)=A(3,1)*A(1,2)-A(3,2)*A(1,1) A(3,3)=A(1,1)*A(2,2)-A(1,2)*A(2,1) RETURN END SUBROUTINE NORM(A) changes made 09-May-1988 *** IMPLICIT NONE C C NORMALIZES COLUMNS OF A C REAL*8 A(3,3),FN INTEGER*4 I,J,K DO 30 I=1,3 FN=0.0D0 DO 10 J=1,3 10 FN=FN+A(J,I)**2 FN=SQRT(FN) DO 20 K=1,3 20 A(K,I)=A(K,I)/FN 30 CONTINUE RETURN END subroutine str_upcase (output,input) c c----convert lower case to upper case implicit none c character*(*) input,output integer length,i c length = len(input) do i=1,length if ((input(i:i) .ge. 'a') .and. (input(i:i) .le. 'z')) then output(i:i) = char(ichar(input(i:i)) - ichar('a') + ichar('A')) endif enddo c c----finish return end subroutine str_trim(out,in,k) character*(*) out,in integer k character*1 tab tab = char(9) k = len(in) out = in i = k do while (i.gt.0 .and. i.eq.k) if (in(k:k).eq.' ' .or. in(k:k).eq.tab) then k = k - 1 end if i = i - 1 end do if (k.gt.len(out)) then k = len(out) end if return end SUBROUTINE LIB$DATE_TIME (DATE_TIME) C C----USE THIS CODE ON A SILICON GRAPHICS MIPS MACHINE IMPLICIT NONE C EXTERNAL TIME INTEGER TIME CHARACTER*24 DATTIM,CTIME CHARACTER*(*) DATE_TIME C DATTIM = CTIME(TIME()) DATE_TIME = DATTIM(9:10)//'-'//DATTIM(5:7)//'-'//DATTIM(21:24)// $ ' '//DATTIM(12:19) C RETURN END