0

I would like to average first column and last column shown here(that is my output). My code does several things such as compare specific atoms with the name "CA" for each frame (I have 1000 frames), exclude the nearest neighbours and depending on the cutoff value, and count those contacts and print them separately according to my wish (input file looks like this). I would like to print a file where it gives me the output something like following as an example:

1 0

2 12

3 12

....

100 16


Need guidance to achieve that by helping me form a loop or a condition.

  open(unit=11,file="0-1000.gro", status="old", action="read")      
  do f=1,frames,10               
  25 format (F10.5,F10.5,F10.5)      

  do h=1,natoms_frames
  read(11,format11)nom(h),resname(h),atmtype(h),num(h),x(h),y(h),z(h)
  end do
  read(11,25)lasta(f),lastb(f),lastc(f) 

  count=0    
  do h=1, natoms_frames        
  if (atmtype(h).eq.'  CA') then
      count=count+1
      CAx(count)=x(h) 
      CAy(count)=y(h)
      CAz(count)=z(h)    
  end if
  end do      

  do h=1, count
  avg_cal=0 
    cal=0 
    do hh=h+3, count
    if (h.ne.hh) then
  ! finding distance formula from the gro file  
  distance = sqrt((CAx(hh)-CAx(h))**2 + (CAy(hh)-CAy(h))**2 + (CAz(hh)-CAz(h))**2)
  if (distance.le.cutoff) then
  cal = cal+1
  set = set+1
  final_set=final_set+1
  avg_cal=avg_cal+1
  end if
  end if
  end do
   write(*,*)h,cal,final_set 
  end do

  end do ! end of frames

  close(11)

  end program num_contacts     
Vishal
  • 3
  • 3
  • 2
    I don't understand your question. That's a lot of code to inspect and I don't think much of it is relevant. I think we need a [mcve] and some sample inputs and outputs -- here on SO, not at the other end of a hyperlink to who-knows-where. – High Performance Mark May 01 '18 at 07:31
  • Please rewrite the double loop in pseudo-code (the first one is ok) and/or better explain the rejection criteria. –  May 02 '18 at 07:15

1 Answers1

0

Writing in a file (in ASCII) is the same as printing on screen, except that you need to specify the file's unit.

So first, open your file:

open(unit = 10, file = "filename.dat", access = "sequential", form = "formatted", status = "new", iostat = ierr)

where 'unit' is a unique integer that will be associated to your file (as an ID). Note that the 'form' is 'formatted' as you want to output something readable by a human. Then you can start your loop:

do h = 1, natoms_frames
  ! BUNCH OF STUFF TO CALCULATE REQUIRED VARIABLES
  write(10, *) h, cal, final_set
end do

Finally, close your file:

close(10)
Kefeng91
  • 802
  • 6
  • 10