add fraction directive

src/prepar/pre_input.F

@@ -1,6 +1,6 @@
       subroutine pre_input(irtdb)
 c
-c $Id: pre_input.F,v 1.40 2000-03-29 18:38:19 d3j191 Exp $
+c $Id: pre_input.F,v 1.41 2000-03-29 20:29:00 d3j191 Exp $
 c
       implicit none
 c
@@ -11,13 +11,19 @@ c
 c
       integer irtdb
       character*255 item,atom,type,atomi,atomj,atomk,atoml
+c
+      integer mcount
+      parameter(mcount=10)
+      integer mfract
+      parameter(mfract=10)
 c
       character*3 slvnam
       character*80 source
       character*10 slvmdl
       character*80 ffield,sysnam
       character*255 dir_s,dir_x,dir_u,dir_t,commnd,modify
-      integer newtop,newseq,newrst,icount,multip
+      integer newtop,newseq,newrst,icount(10),multip,ncount
+      integer nfract,ifract(mfract)
 c
       integer len,numcmd,nummod,mgrid,mnoe,maxscf,num,mset
       integer isgm,icyren,model,nxlnk
@@ -30,7 +36,7 @@ c
       integer nxrep,nyrep,nzrep
       real*8 rrep
       real*8 value1,value2
-      character*4 scount
+      character*4 scount(mcount)
 c
       logical ltouch,lrdcoo,lslvnt,needwt,lwrcoo
 c
@@ -66,7 +72,8 @@ c
       mnoe=0
       maxscf=0
       qscale=-1.0d0
-      icount=0
+      ncount=0
+      nfract=0
       icyren=0
       model=0
       nxlnk=0
@@ -302,13 +309,28 @@ c
       goto 2
       endif
 c
+c     divide molecules into fractions
+c     -------------------------------
+c
+      if(inp_compare(.false.,'fraction',item)) then
+    4 continue
+      if(.not.inp_i(i)) goto 2
+      nfract=nfract+1
+      if(nfract.gt.mfract) call errquit('Too many fractions',nfract)
+      ifract(nfract)=i
+      goto 4
+      endif
+c
 c     add counter ions to sequence file
 c     ---------------------------------
 c
       if(inp_compare(.false.,'counter',item)) then
-      if(.not.inp_i(icount)) call errquit('pre_input: counter',9999)
+      ncount=ncount+1
+      if(ncount.gt.mcount) call errquit('Too many counters',ncount)
+      if(.not.inp_i(icount(ncount)))
+     + call errquit('pre_input: counter',9999)
       if(.not.inp_a(item)) call errquit('pre_input: counter',9999)
-      scount=item(1:4)
+      scount(ncount)=item(1:4)
       goto 2
       endif
 c
@@ -1030,13 +1052,22 @@ c
      + call errquit('pre_input: rtdb_put failed',0)
       if(.not.rtdb_put(irtdb,'prep:newrst',mt_int,1,newrst))
      + call errquit('pre_input: rtdb_put failed',0)
+c
+      if(.not.rtdb_put(irtdb,'prep:nfract',mt_int,1,nfract))
+     + call errquit('pre_input: rtdb_put failed',0)
+      if(nfract.gt.0) then
+      if(.not.rtdb_put(irtdb,'prep:ifract',mt_int,nfract,ifract))
+     + call errquit('pre_input: rtdb_put failed',0)
+      endif
 c
       if(.not.rtdb_put(irtdb,'prep:icyren',mt_int,1,icyren))
      + call errquit('pre_input: rtdb_put failed',0)
-      if(.not.rtdb_put(irtdb,'prep:icount',mt_int,1,icount))
+      if(.not.rtdb_put(irtdb,'prep:ncount',mt_int,1,ncount))
      + call errquit('pre_input: rtdb_put failed',0)
-      if(icount.gt.0) then
-      if(.not.rtdb_cput(irtdb,'prep:scount',1,scount))
+      if(ncount.gt.0) then
+      if(.not.rtdb_put(irtdb,'prep:icount',mt_int,ncount,icount))
+     + call errquit('pre_input: rtdb_put failed',0)
+      if(.not.rtdb_cput(irtdb,'prep:scount',ncount,scount))
      + call errquit('pre_input: rtdb_put failed',0)
       endif
       if(mgrid.gt.0) then

src/prepar/pre_mkseq.F

@@ -2,9 +2,9 @@
      + lfnpdb,filpdb,lfnseq,filseq,lfnpar,filpar,lfnfrg,lfnsgm,
      + lfntmp,filtmp,lfnmod,filmod,dir_s,dir_x,dir_u,dir_t,
      + slvnam,slvmdl,maxscf,qscale,altloc,chain,icyren,model,nxlnk,
-     + icount,scount)
+     + mcount,ncount,icount,scount,mfract,nfract,ifract)
 c
-c $Id: pre_mkseq.F,v 1.32 2000-03-21 23:58:04 d3j191 Exp $
+c $Id: pre_mkseq.F,v 1.33 2000-03-29 20:29:00 d3j191 Exp $
 c
 c     in : integer lfnout    = logical file number output file
 c          char*80 ffield    = force field from [amber]
@@ -36,12 +36,15 @@ c
 c
       integer matm,mseq,mssb,msgm,mbnd,mang,mdih,mimp,mlnk,mato,nato
       integer natm,nseq,nssb,nsgm,nlnk,maxscf,icyren,model,nxlnk
+      integer mcount,ncount
+      integer mfract,nfract,ifract(mfract)
       integer mfrb,i_frb,l_frb
       integer l_lseq,i_lseq,l_cseq,i_cseq
       integer l_lsgm,i_lsgm,l_csgm,i_csgm
       integer l_lssb,i_lssb,l_llnk,i_llnk,l_clnk,i_clnk
       integer l_latm,i_latm,l_catm,i_catm,l_xatm,i_xatm,l_qatm,i_qatm
-      integer l_bnd,i_bnd,l_ang,i_ang,l_dih,i_dih,l_imp,i_imp,icount
+      integer l_bnd,i_bnd,l_ang,i_ang,l_dih,i_dih,l_imp,i_imp
+      integer icount(mcount)
       integer l_lato,i_lato,i_cato,l_cato,l_xato,i_xato,l_qato,i_qato
 c
       real*8 qscale
@@ -54,7 +57,7 @@ c
       character*80 source
       character*10 slvmdl
       character*1 altloc,chain
-      character*4 scount
+      character*4 scount(mcount)
 c
       integer lfnpdb,lfnout,lfnfrg,lfnpar,lfnseq,lfnsgm,lfntmp,lfnmod
       integer igeom,numslv,nlnkf
@@ -383,7 +386,7 @@ c
      + byte_mb(i_cseq),int_mb(i_lseq),mseq,nseq,
      + int_mb(i_lssb),mssb,nssb,
      + int_mb(i_llnk),byte_mb(i_clnk),mlnk,nlnk,nlnkf,slvmdl,
-     + icount,scount))
+     + mcount,ncount,icount,scount,mfract,nfract,ifract))
      + call errquit('pre_wrtseq failed',9999)
 c
       if(util_print('where',print_debug)) then

src/prepar/pre_rtdbin.F

@@ -1,9 +1,11 @@
       logical function pre_rtdbin(irtdb,ffield,
      + dir_s,dir_x,dir_u,dir_t,source,sysnam,calc,slvnam,slvmdl,
-     + newtop,newseq,newrst,icount,mgrid,gdist,mnoe,maxscf,qscale,
-     + altloc,chain,icyren,model,nxlnk,mdold,ignore,scount)
+     + newtop,newseq,newrst,mcount,ncount,icount,mgrid,gdist,mnoe,
+     + maxscf,qscale,
+     + altloc,chain,icyren,model,nxlnk,mdold,ignore,scount,
+     + mfract,nfract,ifract)
 c
-c $Id: pre_rtdbin.F,v 1.18 1999-12-10 23:53:29 d3j191 Exp $
+c $Id: pre_rtdbin.F,v 1.19 2000-03-29 20:29:01 d3j191 Exp $
 c
 c     function to read input for prepare module from rtdb
 c
@@ -15,17 +17,19 @@ c
 #include "rtdb.fh"
 #include "mafdecls.fh"
 c
-      integer irtdb
+      integer irtdb,mcount,ncount
       character*3 slvnam
       character*80 source
       character*80 ffield,sysnam,calc
       character*10 slvmdl
       character*255 dir_s,dir_x,dir_u,dir_t
-      integer newtop,newseq,newrst,icount,icyren,ignore
+      integer newtop,newseq,newrst,icount(mcount),icyren,ignore
       integer len,ndx,mgrid,mnoe,maxscf,model,nxlnk,mdold
+      integer mfract,nfract
+      integer ifract(mfract)
       real*8 gdist,qscale
       character*1 altloc,chain
-      character*4 scount
+      character*4 scount(mcount)
 c
       character*255 key,value
 c
@@ -65,9 +69,18 @@ c
       if(.not.rtdb_get(irtdb,'prep:newseq',mt_int,1,newseq)) newseq=0
       if(.not.rtdb_get(irtdb,'prep:newrst',mt_int,1,newrst)) newrst=0
 c
-      if(.not.rtdb_get(irtdb,'prep:icount',mt_int,1,icount)) icount=0
-      if(icount.gt.0) then
-      if(.not.rtdb_cget(irtdb,'prep:scount',1,scount)) scount='Na  '
+      if(.not.rtdb_get(irtdb,'prep:nfract',mt_int,1,nfract)) nfract=0
+      if(nfract.gt.0) then
+      if(.not.rtdb_get(irtdb,'prep:ifract',mt_int,nfract,ifract))
+     + call errquit('Fraction input problem',0)
+      endif
+c
+      if(.not.rtdb_get(irtdb,'prep:ncount',mt_int,1,ncount)) ncount=0
+      if(ncount.gt.0) then
+      if(.not.rtdb_get(irtdb,'prep:icount',mt_int,ncount,icount))
+     + call errquit('Counter ion input problem',0)
+      if(.not.rtdb_cget(irtdb,'prep:scount',ncount,scount))
+     + call errquit('Counter ion input problem',0)
       endif
       if(.not.rtdb_get(irtdb,'prep:mgrid',mt_int,1,mgrid)) mgrid=24
       if(.not.rtdb_get(irtdb,'prep:rgrid',mt_dbl,1,gdist)) gdist=0.2d0

src/prepar/pre_wrtseq.F

@@ -1,8 +1,8 @@
       logical function pre_wrtseq(lfnseq,filseq,
      + cseq,lseq,mseq,nseq,lssb,mssb,nssb,llnk,clnk,mlnk,nlnk,nlnkf,
-     + slvmdl,icount,scount)
+     + slvmdl,mcount,ncount,icount,scount,mfract,nfract,ifract)
 c
-c $Id: pre_wrtseq.F,v 1.9 2000-02-20 16:05:15 d3j191 Exp $
+c $Id: pre_wrtseq.F,v 1.10 2000-03-29 20:29:01 d3j191 Exp $
 c
 c     function to write the sequence file
 c
@@ -14,8 +14,10 @@ c
       character*5 cseq(2,mseq)
       character*10 slvmdl
       character*255 filseq
-      integer i,length,mol,icount,ll
-      character*4 scount
+      integer mcount,ncount
+      integer i,j,l,length,mol,molec,icount(mcount),ll
+      character*4 scount(mcount)
+      integer mfract,nfract,ifract(mfract)
 c
       length=index(filseq,' ')-1
 c
@@ -24,10 +26,23 @@ c
 c
       ll=0
       mol=0
+      molec=1
       do 1 i=1,nseq
       if(i.gt.1.and.lseq(4,i).ne.mol) then
+      if(nfract.eq.0) then
       write(lfnseq,1001)
  1001 format('molecule')
+      else
+      do 11 j=1,nfract
+      if(ifract(j).eq.molec) then
+      write(lfnseq,1005)
+      goto 111
+      endif
+   11 continue
+      write(lfnseq,1001)
+  111 continue
+      endif
+      molec=molec+1
       endif
       mol=lseq(4,i)
       length=index(cseq(2,i),'_')-1
@@ -38,14 +53,18 @@ c      write(lfnseq,1000) lseq(1,i),cseq(2,i)(1:length)
       ll=lseq(1,i)
     1 continue
 c
-      if(icount.gt.0) then
+      if(ncount.gt.0) then
       write(lfnseq,1005)
  1005 format('fraction')
-      do 3 i=1,icount
-      write(lfnseq,1006) ll+i,scount
+      l=0
+      do 3 i=1,ncount
+      do 33 j=1,icount(i)
+      l=l+1
+      write(lfnseq,1006) ll+l,scount(i)
  1006 format(i5,a)
-      if(i.lt.icount) write(lfnseq,1007)
+      if(i.lt.ncount.or.j.lt.icount(i)) write(lfnseq,1007)
  1007 format('molecule')
+   33 continue
     3 continue
       endif
 c

src/prepar/prepar.F

@@ -1,6 +1,6 @@
        logical function prepar(irtdb0)
 c
-c $Id: prepar.F,v 1.32 1999-11-04 22:53:08 d3j191 Exp $
+c $Id: prepar.F,v 1.33 2000-03-29 20:29:01 d3j191 Exp $
 c
 c ********************************************************
 c ********************************************************
@@ -39,6 +39,11 @@ c
       external pre_mknoe
 c
       integer irtdb0,irtdb,itask
+c
+      integer mcount
+      parameter(mcount=10)
+      integer mfract
+      parameter(mfract=10)
 c
       character*255 filpdb,filseq,filtop,filrst,filpar,filtmp,filcmd
       character*255 filmod,prefix,filxyz,filqqq,filnoe,filpmf
@@ -48,12 +53,14 @@ c
       character*3 slvnam
       character*80 source
       character*1 altloc,chain
-      character*4 scount
+      character*4 scount(mcount)
 c
       integer lfnpdb,lfnout,lfnfrg,lfnseq,lfnsgm,lfntop,lfnrst,lfnpar
       integer lfnxyz,lfnqqq,lfnnoe,lfnpmf
       integer len,lenc,lend,lfntmp,lfncmd,lfnslv,lfnmod,nxlnk,mdold
-      integer newtop,newseq,newrst,mgrid,mnoe,maxscf,icount,icyren,model
+      integer newtop,newseq,newrst,mgrid,mnoe,maxscf,icount(mcount)
+      integer ncount,icyren,model
+      integer nfract,ifract(mfract)
       integer ignore
       real*8 gdist,qscale
       logical lstate
@@ -112,9 +119,10 @@ c     get info from rtdb
 c     ------------------
 c
       if(.not.pre_rtdbin(irtdb,ffield,dir_s,dir_x,dir_u,dir_t,source,
-     + sysnam,calc,slvnam,slvmdl,newtop,newseq,newrst,icount,mgrid,
+     + sysnam,calc,slvnam,slvmdl,newtop,newseq,newrst,mcount,ncount,
+     + icount,mgrid,
      + gdist,mnoe,maxscf,qscale,altloc,chain,icyren,model,nxlnk,mdold,
-     + ignore,scount))
+     + ignore,scount,mfract,nfract,ifract))
      + call errquit('pre_rtdbin failed',9999)
 c
 c     directories
@@ -177,7 +185,6 @@ c
       if(source(1:1).eq.' ') source='geometry'
       if(source(1:4).eq.'rtdb') source='geometry'
 c
-c
 c     check if the topology file exists
 c     ---------------------------------
 c
@@ -229,7 +236,7 @@ c
      + lfnpdb,filpdb,lfnseq,filseq,lfnpar,filpar,lfnfrg,lfnsgm,
      + lfntmp,filtmp,lfnmod,filmod,dir_s,dir_x,dir_u,dir_t,
      + slvnam,slvmdl,maxscf,qscale,altloc,chain,icyren,model,nxlnk,
-     + icount,scount))
+     + mcount,ncount,icount,scount,mfract,nfract,ifract))
      + call errquit('pre_mkseq failed',9999)
       if(util_print('files',print_default)) then
       write(lfnout,2009) filseq(1:index(filseq,' ')-1)