Changeset bd2278d for EXAMPLES


Ignore:
Timestamp:
09/05/08 11:49:42 (16 years ago)
Author:
baerbaer <baerbaer@…>
Branches:
master
Children:
fafe4d6
Parents:
2ebb8b6
Message:

Reformatting comments and continuation marks.

Fortran 90 and higher use ! to mark comments no matter where they are in the
code. The only valid continuation marker is &.
I also added the SMMP.kdevelop.filelist to the repository to make it easier
to use kdevelop.

git-svn-id: ​svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@12 26dc1dd8-5c4e-0410-9ffe-d298b4865968

Location:
EXAMPLES
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • EXAMPLES/parallel_tempering_s.f

    r2ebb8b6 rbd2278d  
    2020      integer switch
    2121
    22 c =================================================== Energy setup
     22! =================================================== Energy setup
    2323
    24 c            Directory for SMMP libraries
    25 c     Change the following directory path to where you want to put SMMP
    26 c     libraries of residues.
     24!            Directory for SMMP libraries
     25!     Change the following directory path to where you want to put SMMP
     26!     libraries of residues.
    2727      libdir='../SMMP/'
    2828
    29 c      The switch in the following line is now not used.
     29!      The switch in the following line is now not used.
    3030      flex=.false.        ! .true. for Flex  / .false. for ECEPP
    3131
    32 c     Choose energy type with the following switch instead ...
     32!     Choose energy type with the following switch instead ...
    3333      ientyp = 0
    34 c        0  => ECEPP2 or ECEPP3 depending on the value of sh2
    35 c        1  => FLEX
    36 c        2  => Lund force field
    37 c        3  => ECEPP with Abagyan corrections
    38 c
     34!        0  => ECEPP2 or ECEPP3 depending on the value of sh2
     35!        1  => FLEX
     36!        2  => Lund force field
     37!        3  => ECEPP with Abagyan corrections
     38!
    3939
    4040      sh2=.false.         ! .true. for ECEPP/2; .false. for ECEPP3
    … …  
    4848      call init_energy(libdir)
    4949
    50 c ================================================= Structure setup
     50! ================================================= Structure setup
    5151
    5252      grpn = 'nh2' ! N-terminal group
    … …  
    5959      ntlml = 0
    6060      write (*,*) 'Solvent: ', itysol
    61 c     Initialize random number generator.
     61!     Initialize random number generator.
    6262      call sgrnd(31433)
    6363     
    … …  
    6969
    7070      call init_molecule(iabin,grpn,grpc,seqfile,varfile)
    71 c Decide if and when to use BGS, and initialize Lund data structures
     71! Decide if and when to use BGS, and initialize Lund data structures
    7272      bgsprob=0.75   ! Prob for BGS, given that it is possible
    73 c upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
    74 c 2 => temperature dependent choice
     73! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
     74! 2 => temperature dependent choice
    7575      upchswitch=1
    7676      rndord=.true.
    … …  
    8080     
    8181
    82 c ========================================  Add your task down here
     82! ========================================  Add your task down here
    8383      num_rep = 5
    8484      nequi = 100
    … …  
    8787      newsta = .true.
    8888      switch = 1
    89 c     parallel tempering on a single CPU
     89!     parallel tempering on a single CPU
    9090      eol = energy()
    9191      write (*,*) "Energy before randomization:", eol
    … …  
    9494      write (*,*) "Final energy:", eol
    9595
    96 c ========================================  End of main     
     96! ========================================  End of main     
    9797       end

  • EXAMPLES/partem_p.f

    r2ebb8b6 rbd2278d  
    1 c**************************************************************
    2 c     
    3 c This file contains the subroutines: partem_p
    4 C Compared to the version in the main distribution, this
    5 C routine doesn't write the rmsd nor native contacts to the time
    6 C series.
    7 c     
    8 c Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    9 c                      Shura Hayryan, Chin-Ku Hu
    10 c Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    11 c                      Jan H. Meinke, Sandipan Mohanty
    12 c     
    13 c     **************************************************************
     1!**************************************************************
     2!     
     3! This file contains the subroutines: partem_p
     4! Compared to the version in the main distribution, this
     5! routine doesn't write the rmsd nor native contacts to the time
     6! series.
     7!     
     8! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
     9!                      Shura Hayryan, Chin-Ku Hu
     10! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
     11!                      Jan H. Meinke, Sandipan Mohanty
     12!     
     13!     **************************************************************
    1414
    1515      subroutine  partem_p(num_rep, nequi, nswp, nmes, nsave, newsta,
    1616     &                     switch, rep_id, partem_comm)
    17 C     
    18 C     PURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM
    19 C     ON PARALLEL COMPUTERS USING MPI
    20 C     
    21 C     switch: Choses the starting configuration:
    22 C     -1 - stretched configuration
    23 C     0 - don't change anything
    24 C     1 - random start configuration
    25 C     
    26 c     CALLS:  addang,contacts,energy,hbond,helix,iendst,metropolis,
    27 c     outvar,(rand),rgyr
    28 C     
     17!     
     18!     PURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM
     19!     ON PARALLEL COMPUTERS USING MPI
     20!     
     21!     switch: Choses the starting configuration:
     22!     -1 - stretched configuration
     23!     0 - don't change anything
     24!     1 - random start configuration
     25!     
     26!     CALLS:  addang,contacts,energy,hbond,helix,iendst,metropolis,
     27!     outvar,(rand),rgyr
     28!     
    2929      include '../INCL.H'
    3030      include '../INCP.H'
    … …  
    3737      external can_weight
    3838
    39 C     nequi:  number of Monte Carlo sweeps for thermalization
    40 C     nswp:   number of Monte Carlo sweeps
    41 C     nmes:   number of Monte Carlo sweeps between measurments
    42 C     newsta: .true. for new simulations, .false. for re-start
     39!     nequi:  number of Monte Carlo sweeps for thermalization
     40!     nswp:   number of Monte Carlo sweeps
     41!     nmes:   number of Monte Carlo sweeps between measurments
     42!     newsta: .true. for new simulations, .false. for re-start
    4343
    4444      dimension  eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC),
    … …  
    5050      double precision    e_min, e_minp(MAX_PROC), e_minpt(MAX_PROC)
    5151      integer   h_max, h_maxp(MAX_PROC)
    52 c     Order of replica exchange
     52!     Order of replica exchange
    5353      integer   odd
    5454!     Counter to keep random number generators in sync
    5555      integer randomCount
    5656     
    57 c     Collect partial energies. Only the root writes to disk. We have to
    58 c     collect the information from the different replicas and provide
    59 c     arrays to store them.
    60 c     eyslr    storage array for solvent energy
    61 c     eyelp     -      "        - coulomb energy
    62 c     eyvwp     -      "        - van-der-Waals energy
    63 c     eyhbp     -      "        - hydrogen bonding energy
    64 c     eysmi    -      "        - intermolecular interaction energy
     57!     Collect partial energies. Only the root writes to disk. We have to
     58!     collect the information from the different replicas and provide
     59!     arrays to store them.
     60!     eyslr    storage array for solvent energy
     61!     eyelp     -      "        - coulomb energy
     62!     eyvwp     -      "        - van-der-Waals energy
     63!     eyhbp     -      "        - hydrogen bonding energy
     64!     eysmi    -      "        - intermolecular interaction energy
    6565      double precision eyslr(MAX_PROC)
    6666      double precision eyelp(MAX_PROC),eyvwp(MAX_PROC),eyhbp(MAX_PROC),
    6767     &     eyvrp(MAX_PROC),eysmip(MAX_PROC)
    68 c     Collect information about accessible surface and van-der-Waals volume
    69 c     asap      storage array for solvent accessible surface
    70 c     vdvolp     storage array for van-der-Waals volume
     68!     Collect information about accessible surface and van-der-Waals volume
     69!     asap      storage array for solvent accessible surface
     70!     vdvolp     storage array for van-der-Waals volume
    7171      double precision asap(MAX_PROC), vdvolp(MAX_PROC)
    7272
    … …  
    7575      integer imhbp(MAX_PROC)
    7676      character*80 filebase, fileNameMP, tbase0,tbase1
    77 c     frame     frame number for writing configurations
    78 c     trackID   configuration that should be tracked and written out
    79 c     dir          direction in random walk
    80 c     -1 - visited highest temperature last
    81 c     1 - visited lowest temperature last
    82 c     0 - haven't visited the boundaries yet.
    83 c     dirp      storage array for directions.
     77!     frame     frame number for writing configurations
     78!     trackID   configuration that should be tracked and written out
     79!     dir          direction in random walk
     80!     -1 - visited highest temperature last
     81!     1 - visited lowest temperature last
     82!     0 - haven't visited the boundaries yet.
     83!     dirp      storage array for directions.
    8484      integer frame, trackID, dir
    8585      integer dirp(MAX_PROC)
    … …  
    9292     &            rep_id, num_rep, partem_comm, myrank
    9393      call flush(6)
    94 C     
    95 c     
    96 C     File with temperatures
     94!     
     95!     
     96!     File with temperatures
    9797      open(11,file='temperatures_abeta',status='old')
    9898     
    … …  
    100100      open(18,file=fileNameMP(tbase0,5,9,rep_id),status='unknown')
    101101      if (rep_id.eq.0.and.myrank.eq.0) then
    102 c     File with time series of simulation
     102!     File with time series of simulation
    103103         open(14,file='ts.d',status='unknown')
    104104      endif
    105105     
    106 C     READ IN TEMPERATURES
     106!     READ IN TEMPERATURES
    107107      do i=1,num_rep
    108108         read(11,*) j,temp
    … …  
    111111      close(11)
    112112
    113 c     nresi:  number of residues
     113!     nresi:  number of residues
    114114      nresi=irsml2(1)-irsml1(1)+1
    115 C     
    116 C     Initialize variables
     115!     
     116!     Initialize variables
    117117      do i=1,num_rep     
    118118         acx1(i) = 0.0d0
    … …  
    132132      dir = dirp(rep_id + 1)
    133133
    134 c     _________________________________ Initialize Variables
     134!     _________________________________ Initialize Variables
    135135      if(newsta) then
    136136         iold=0
    … …  
    139139            intem(i) = i
    140140         end do
    141 c     _________________________________ initialize starting configuration
     141!     _________________________________ initialize starting configuration
    142142         if (switch.ne.0) then
    143143            do i=1,nvr
    … …  
    173173         CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
    174174         CALL MPI_BCAST(YOL,num_rep,MPI_DOUBLE_PRECISION,0,
    175      #        MPI_COMM_WORLD,IERR)
     175     &        MPI_COMM_WORLD,IERR)
    176176         CALL MPI_BCAST(E_MINP, num_rep, MPI_DOUBLE_PRECISION, 0,
    177      #        MPI_COMM_WORLD, IERR)
     177     &        MPI_COMM_WORLD, IERR)
    178178         CALL MPI_BCAST(h_maxp,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
    179      $        IERR)
     179     &        IERR)
    180180      end if
    181181     
    … …  
    189189         write(*,*) rep_id, yol(rep_id + 1), eol
    190190      endif
    191 C     Start of simulation
     191!     Start of simulation
    192192      write (*,*) '[',rep_id, myrank, beta, partem_comm,
    193193     &            '] Energy before equilibration:', eol
    194 c     =====================Equilibration by canonical Metropolis
     194!     =====================Equilibration by canonical Metropolis
    195195      do nsw=1,nequi
    196196         call metropolis(eol,acz,can_weight)
    … …  
    199199      write (*,*) '[',rep_id,'] Energy after equilibration:', eol
    200200      call flush(6)
    201 C     
    202 C======================Multiple Markov Chains
     201!     
     202!======================Multiple Markov Chains
    203203      acz = 0
    204204      do nsw=1,nswp
    205 c------------First ordinary Metropolis
     205!------------First ordinary Metropolis
    206206         call metropolis(eol,acz,can_weight)
    207207         iold = iold + 1       
    … …  
    214214            endif
    215215            acz0 = acz
    216 C     Measure global radius of gyration
     216!     Measure global radius of gyration
    217217            call rgyr(0,rgy,ee) 
    218218            rgyp = rgy
    219 C     Measure Helicity and Sheetness
     219!     Measure Helicity and Sheetness
    220220            call helix(nhel,mhel,nbet,mbet)
    221 C     Measure Number of hydrogen bonds
     221!     Measure Number of hydrogen bonds
    222222            mhb = 0
    223223            do i = 1, ntlml
    … …  
    226226            enddo
    227227            call interhbond(imhb)
    228 C     Measure total number of contacts (NCTOT) and number of
    229 C     native contacts (NCNAT)
     228!     Measure total number of contacts (NCTOT) and number of
     229!     native contacts (NCNAT)
    230230            call contacts(nctot,ncnat,dham)
    231 c     Add tracking of lowest energy configuration
     231!     Add tracking of lowest energy configuration
    232232            if (eol.lt.e_min) then
    233 c     Write out configuration
     233!     Write out configuration
    234234               i=rep_id+1
    235235               j=inode(i)
    … …  
    248248               close(15)
    249249            endif
    250 c     Add tracking of configuration with larges hydrogen contents.
     250!     Add tracking of configuration with larges hydrogen contents.
    251251            if ((mhb + imhb).gt.h_max) then
    252 c     Write out configuration
     252!     Write out configuration
    253253               i = rep_id + 1
    254254               j = inode(i)
    … …  
    268268            endif
    269269
    270 C     
    271 C--------------------Gather measurement data
     270!     
     271!--------------------Gather measurement data
    272272! I only use the master node of each replica for data collection. The
    273273! variable partem_comm provides the appropriate communicator.
    … …  
    310310     &              MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
    311311
    312 c     Write trajectory
     312!     Write trajectory
    313313               write (18,*) '@@@',iold,inode(rep_id+1)
    314314               call outvbs(0,18)
    315315               write (18,*) '###'
    316316!                call flush(18)
    317 c     Write current configuration
     317!     Write current configuration
    318318               if ((mod(iold, nsave).eq.0)) then
    319319                  filebase = "conf_0000.var"
    … …  
    324324            if(rep_id.eq.0.and.myrank.eq.0) then
    325325               randomCount = 0
    326 c  Update acceptance, temperature wise average of E and E^2 used to calculate
    327 c  specific heat.
     326!  Update acceptance, temperature wise average of E and E^2 used to calculate
     327!  specific heat.
    328328               do i=1,num_rep
    329329                  j=intem(i)
    330330                  acy(i)=0.0
    331 c  Above: contents of acy1 are added to acy(i) a few lines down.
    332 c  acy1(intem(i)) contains information received from the node at temperature
    333 c  i, on how many updates have been accepted in node intem(i). Since acz
    334 c  is not reset to 0 every cycle, acy(i) must be set to 0 here. Else, there
    335 c  will be serious double counting and the values of acceptance printed
    336 c  will be simply wrong.
     331!  Above: contents of acy1 are added to acy(i) a few lines down.
     332!  acy1(intem(i)) contains information received from the node at temperature
     333!  i, on how many updates have been accepted in node intem(i). Since acz
     334!  is not reset to 0 every cycle, acy(i) must be set to 0 here. Else, there
     335!  will be serious double counting and the values of acceptance printed
     336!  will be simply wrong.
    337337                  e_minpt(i)=e_minp(intem(i))
    338338               end do
    … …  
    346346
    347347
    348 C     Write measurements to the time series file ts.d
     348!     Write measurements to the time series file ts.d
    349349               do i=1,num_rep
    350350                  j=intem(i)
    … …  
    354354                     
    355355               end do
    356 c     Write the current parallel tempering information into par_R.in
     356!     Write the current parallel tempering information into par_R.in
    357357               if ((mod(iold, nsave).eq.0))
    358358     &         then
    … …  
    363363     &                    h_maxp(i)
    364364                  end do
    365 C     -------------------------- Various statistics of current run
    366 c               swp=nswp-nequi
     365!     -------------------------- Various statistics of current run
     366!               swp=nswp-nequi
    367367                  swp=nsw
    368368                  write(13,*) 'Acceptance rate for change of chains:'
    … …  
    370370                     temp=1.0d0/pbe(k1)/0.00198773
    371371                     write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
    372 c  Above: it's the acceptance rate of exchange of replicas. Since a
    373 c  replica exchange is attempted only once every nmes sweeps, the
    374 c  rate should be normalized with (nmes/swp).
     372!  Above: it's the acceptance rate of exchange of replicas. Since a
     373!  replica exchange is attempted only once every nmes sweeps, the
     374!  rate should be normalized with (nmes/swp).
    375375                  end do
    376376                  write(13,*)
    … …  
    381381                     geavm(k1) = nmes*eavm(k1)/swp
    382382                     gsph(k1)  = (nmes*sph(k1)/swp-geavm(k1)**2)
    383      #                    *beta*beta/nresi
     383     &                    *beta*beta/nresi
    384384                     write(13,'(a,2f9.2,i4,f12.3)')
    385385     &                    'Temperature, Node,local acceptance rate:',
    386386     &                    beta,temp,k,acy(k1)/dble(nsw*nvr)
    387 c  Above: Changed (nswp-nequi) in the denominator of acceptance as
    388 c  acceptance values are initialized to 0 after equilibration cycles are
    389 c  finished. Note also that since this is being written in the middle of
    390 c  the simulation, it is normalized to nsw instead of nswp.
     387!  Above: Changed (nswp-nequi) in the denominator of acceptance as
     388!  acceptance values are initialized to 0 after equilibration cycles are
     389!  finished. Note also that since this is being written in the middle of
     390!  the simulation, it is normalized to nsw instead of nswp.
    391391                     write(13,'(a,3f12.2)')
    392392     &                    'Last Energy, Average Energy, Spec. Heat:',
    … …  
    401401               end if   
    402402
    403 C--------------------Parallel Tempering  update
    404 c     Swap with right neighbor (odd, even)           
     403!--------------------Parallel Tempering  update
     404!     Swap with right neighbor (odd, even)           
    405405               if(odd.eq.1) then
    406406                  nu=1
    407407                  no1 = num_rep-1
    408 c     Swap with left neighbor (even, odd)
     408!     Swap with left neighbor (even, odd)
    409409               else
    410410                  nu = 2
    … …  
    413413               do i=nu,no1,2
    414414                  j=i+1
    415 c     Periodic bc for swaps
     415!     Periodic bc for swaps
    416416                  if(i.eq.num_rep) j=1
    417417                  iidx=intem(i)
    … …  
    429429                  end if
    430430               end do
    431 c     ---------------- End Loop over nodes which creates a new temperature
    432 c     map for all nodes, at the node with rank 0.
    433 c     
     431!     ---------------- End Loop over nodes which creates a new temperature
     432!     map for all nodes, at the node with rank 0.
     433!     
    434434               odd = 1 - odd
    435435            end if
    436 c     End of "if (myrank.eq.0) ...". The block above includes PT update and
    437 c     writing of observables into the time series file etc.
     436!     End of "if (myrank.eq.0) ...". The block above includes PT update and
     437!     writing of observables into the time series file etc.
    438438           
    439 c     Below: Communicate new temperature-node map to all nodes
     439!     Below: Communicate new temperature-node map to all nodes
    440440            CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
    441441     &           IERR)
    … …  
    446446            CALL MPI_BCAST(H_MAXP,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
    447447     &           IERR)
    448 c Synchronize random number generators for replica 0
     448! Synchronize random number generators for replica 0
    449449            if (rep_id.eq.0) then
    450450               CALL MPI_BCAST(randomCount,1,MPI_INTEGER,0,my_mpi_comm,
    … …  
    467467
    468468         endif
    469 c        End of "if (mod(iold,nmes).eq.0) ..."
     469!        End of "if (mod(iold,nmes).eq.0) ..."
    470470      end do
    471 c-----------End Loop over sweeps
    472 c     
    473 C     OUTPUT:
    474 C--------------------For Re-starts:
     471!-----------End Loop over sweeps
     472!     
     473!     OUTPUT:
     474!--------------------For Re-starts:
    475475      nu = rep_id + 1
    476476      filebase = "conf_0000.var"
    … …  
    484484      if (partem_comm.ne.MPI_COMM_NULL) then
    485485         CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,
    486      #        MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
     486     &        MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
    487487         CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1,
    488      #     MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
     488     &     MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
    489489      endif
    490490     
    … …  
    496496            write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
    497497         end do
    498 C     -------------------------- Various statistics of current run
     498!     -------------------------- Various statistics of current run
    499499         swp=nswp
    500500         write(13,*) 'Acceptance rate for change of chains:'
    … …  
    519519         end do
    520520         close(13)
    521 c         close(16)
     521!         close(16)
    522522      end if
    523523      close(18)
    524524
    525 c     =====================
     525!     =====================
    526526      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
    527527
Note: See TracChangeset for help on using the changeset viewer.