8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 integer :: ichir1,ichir2
748 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
749 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
750 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 ! For printing parameters after they are read set the following in the UNRES
755 ! setenv PRINT_PARM YES
757 ! To print parameters in LaTeX format rather than as ASCII tables:
761 call getenv_loc("PRINT_PARM",lancuch)
762 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
763 call getenv_loc("LATEX",lancuch)
764 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
766 dwa16=2.0d0**(1.0d0/6.0d0)
768 ! Assign virtual-bond length
771 vblinv2=vblinv*vblinv
773 ! Read the virtual-bond parameters, masses, and moments of inertia
774 ! and Stokes' radii of the peptide group and side chains
776 allocate(dsc(ntyp1)) !(ntyp1)
777 allocate(dsc_inv(ntyp1)) !(ntyp1)
778 allocate(nbondterm(ntyp)) !(ntyp)
779 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
780 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(msc(ntyp+1)) !(ntyp+1)
782 allocate(isc(ntyp+1)) !(ntyp+1)
783 allocate(restok(ntyp+1)) !(ntyp+1)
784 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
790 read (ibond,*) vbldp0,akp,mp,ip,pstok
793 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
794 dsc(i) = vbldsc0(1,i)
798 dsc_inv(i)=1.0D0/dsc(i)
802 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
804 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
805 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
806 dsc(i) = vbldsc0(1,i)
810 dsc_inv(i)=1.0D0/dsc(i)
815 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
816 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
818 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
820 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
821 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
823 write (iout,'(13x,3f10.5)') &
824 vbldsc0(j,i),aksc(j,i),abond0(j,i)
828 !----------------------------------------------------
829 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
830 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
831 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
832 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
833 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
834 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
847 ! Read the parameters of the probability distribution/energy expression
848 ! of the virtual-bond valence angles theta
851 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
852 (bthet(j,i,1,1),j=1,2)
853 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
854 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
855 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
859 athet(1,i,1,-1)=athet(1,i,1,1)
860 athet(2,i,1,-1)=athet(2,i,1,1)
861 bthet(1,i,1,-1)=-bthet(1,i,1,1)
862 bthet(2,i,1,-1)=-bthet(2,i,1,1)
863 athet(1,i,-1,1)=-athet(1,i,1,1)
864 athet(2,i,-1,1)=-athet(2,i,1,1)
865 bthet(1,i,-1,1)=bthet(1,i,1,1)
866 bthet(2,i,-1,1)=bthet(2,i,1,1)
870 athet(1,i,-1,-1)=athet(1,-i,1,1)
871 athet(2,i,-1,-1)=-athet(2,-i,1,1)
872 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
873 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
874 athet(1,i,-1,1)=athet(1,-i,1,1)
875 athet(2,i,-1,1)=-athet(2,-i,1,1)
876 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
877 bthet(2,i,-1,1)=bthet(2,-i,1,1)
878 athet(1,i,1,-1)=-athet(1,-i,1,1)
879 athet(2,i,1,-1)=athet(2,-i,1,1)
880 bthet(1,i,1,-1)=bthet(1,-i,1,1)
881 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
886 polthet(j,i)=polthet(j,-i)
889 gthet(j,i)=gthet(j,-i)
897 'Parameters of the virtual-bond valence angles:'
898 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
899 ' ATHETA0 ',' A1 ',' A2 ',&
902 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
903 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
905 write (iout,'(/a/9x,5a/79(1h-))') &
906 'Parameters of the expression for sigma(theta_c):',&
907 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
908 ' ALPH3 ',' SIGMA0C '
910 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
911 (polthet(j,i),j=0,3),sigc0(i)
913 write (iout,'(/a/9x,5a/79(1h-))') &
914 'Parameters of the second gaussian:',&
915 ' THETA0 ',' SIGMA0 ',' G1 ',&
918 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
919 sig0(i),(gthet(j,i),j=1,3)
923 'Parameters of the virtual-bond valence angles:'
924 write (iout,'(/a/9x,5a/79(1h-))') &
925 'Coefficients of expansion',&
926 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
927 ' b1*10^1 ',' b2*10^1 '
929 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
930 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
931 (10*bthet(j,i,1,1),j=1,2)
933 write (iout,'(/a/9x,5a/79(1h-))') &
934 'Parameters of the expression for sigma(theta_c):',&
935 ' alpha0 ',' alph1 ',' alph2 ',&
936 ' alhp3 ',' sigma0c '
938 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
939 (polthet(j,i),j=0,3),sigc0(i)
941 write (iout,'(/a/9x,5a/79(1h-))') &
942 'Parameters of the second gaussian:',&
943 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
946 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
947 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
953 ! Read the parameters of Utheta determined from ab initio surfaces
954 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
956 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
957 ntheterm3,nsingle,ndouble
958 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
960 !----------------------------------------------------
961 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
962 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
963 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
964 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
965 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
966 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
967 !(maxtheterm,-maxthetyp1:maxthetyp1,&
968 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
969 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
970 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
971 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
972 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
973 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
974 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
975 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
976 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
977 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
978 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
979 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
980 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
981 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
982 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
983 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
984 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
986 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
988 ithetyp(i)=-ithetyp(-i)
991 aa0thet(:,:,:,:)=0.0d0
992 aathet(:,:,:,:,:)=0.0d0
993 bbthet(:,:,:,:,:,:)=0.0d0
994 ccthet(:,:,:,:,:,:)=0.0d0
995 ddthet(:,:,:,:,:,:)=0.0d0
996 eethet(:,:,:,:,:,:)=0.0d0
997 ffthet(:,:,:,:,:,:,:)=0.0d0
998 ggthet(:,:,:,:,:,:,:)=0.0d0
1000 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1002 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1003 ! VAR:1=non-glicyne non-proline 2=proline
1004 ! VAR:negative values for D-aminoacid
1006 do j=-nthetyp,nthetyp
1007 do k=-nthetyp,nthetyp
1008 read (ithep,'(6a)',end=111,err=111) res1
1009 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1010 ! VAR: aa0thet is variable describing the average value of Foureir
1011 ! VAR: expansion series
1012 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1013 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1014 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1015 read (ithep,*,end=111,err=111) &
1016 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1017 read (ithep,*,end=111,err=111) &
1018 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1019 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1020 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1021 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1023 read (ithep,*,end=111,err=111) &
1024 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1025 ffthet(lll,llll,ll,i,j,k,iblock),&
1026 ggthet(llll,lll,ll,i,j,k,iblock),&
1027 ggthet(lll,llll,ll,i,j,k,iblock),&
1028 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1033 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1034 ! coefficients of theta-and-gamma-dependent terms are zero.
1035 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1036 ! RECOMENTDED AFTER VERSION 3.3)
1040 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1041 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1043 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1044 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1047 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1049 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1052 ! AND COMMENT THE LOOPS BELOW
1056 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1057 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1059 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1060 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1063 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1065 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1070 ! Substitution for D aminoacids from symmetry.
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1077 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1081 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1082 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1083 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1084 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1090 ffthet(llll,lll,ll,i,j,k,iblock)= &
1091 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1092 ffthet(lll,llll,ll,i,j,k,iblock)= &
1093 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1094 ggthet(llll,lll,ll,i,j,k,iblock)= &
1095 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1096 ggthet(lll,llll,ll,i,j,k,iblock)= &
1097 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1106 ! Control printout of the coefficients of virtual-bond-angle potentials
1109 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1114 write (iout,'(//4a)') &
1115 'Type ',onelett(i),onelett(j),onelett(k)
1116 write (iout,'(//a,10x,a)') " l","a[l]"
1117 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1118 write (iout,'(i2,1pe15.5)') &
1119 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1121 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1122 "b",l,"c",l,"d",l,"e",l
1124 write (iout,'(i2,4(1pe15.5))') m,&
1125 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1126 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1130 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1131 "f+",l,"f-",l,"g+",l,"g-",l
1134 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1135 ffthet(n,m,l,i,j,k,iblock),&
1136 ffthet(m,n,l,i,j,k,iblock),&
1137 ggthet(n,m,l,i,j,k,iblock),&
1138 ggthet(m,n,l,i,j,k,iblock)
1148 write (2,*) "Start reading THETA_PDB",ithep_pdb
1150 ! write (2,*) 'i=',i
1151 read (ithep_pdb,*,err=111,end=111) &
1152 a0thet(i),(athet(j,i,1,1),j=1,2),&
1153 (bthet(j,i,1,1),j=1,2)
1154 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1155 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1156 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1157 sigc0(i)=sigc0(i)**2
1160 athet(1,i,1,-1)=athet(1,i,1,1)
1161 athet(2,i,1,-1)=athet(2,i,1,1)
1162 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1163 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1164 athet(1,i,-1,1)=-athet(1,i,1,1)
1165 athet(2,i,-1,1)=-athet(2,i,1,1)
1166 bthet(1,i,-1,1)=bthet(1,i,1,1)
1167 bthet(2,i,-1,1)=bthet(2,i,1,1)
1170 a0thet(i)=a0thet(-i)
1171 athet(1,i,-1,-1)=athet(1,-i,1,1)
1172 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1173 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1174 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1175 athet(1,i,-1,1)=athet(1,-i,1,1)
1176 athet(2,i,-1,1)=-athet(2,-i,1,1)
1177 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1178 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1179 athet(1,i,1,-1)=-athet(1,-i,1,1)
1180 athet(2,i,1,-1)=athet(2,-i,1,1)
1181 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1182 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1183 theta0(i)=theta0(-i)
1187 polthet(j,i)=polthet(j,-i)
1190 gthet(j,i)=gthet(j,-i)
1193 write (2,*) "End reading THETA_PDB"
1198 !-------------------------------------------
1199 allocate(nlob(ntyp1)) !(ntyp1)
1200 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1201 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1202 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1209 gaussc(:,:,:,:)=0.0D0
1213 ! Read the parameters of the probability distribution/energy expression
1214 ! of the side chains.
1217 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1221 dsc_inv(i)=1.0D0/dsc(i)
1232 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1233 ((blower(k,l,1),l=1,k),k=1,3)
1234 censc(1,1,-i)=censc(1,1,i)
1235 censc(2,1,-i)=censc(2,1,i)
1236 censc(3,1,-i)=-censc(3,1,i)
1238 read (irotam,*,end=112,err=112) bsc(j,i)
1239 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1240 ((blower(k,l,j),l=1,k),k=1,3)
1241 censc(1,j,-i)=censc(1,j,i)
1242 censc(2,j,-i)=censc(2,j,i)
1243 censc(3,j,-i)=-censc(3,j,i)
1244 ! BSC is amplitude of Gaussian
1251 akl=akl+blower(k,m,j)*blower(l,m,j)
1255 if (((k.eq.3).and.(l.ne.3)) &
1256 .or.((l.eq.3).and.(k.ne.3))) then
1257 gaussc(k,l,j,-i)=-akl
1258 gaussc(l,k,j,-i)=-akl
1260 gaussc(k,l,j,-i)=akl
1261 gaussc(l,k,j,-i)=akl
1270 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1273 if (nlobi.gt.0) then
1275 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1276 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1277 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1278 'log h',(bsc(j,i),j=1,nlobi)
1279 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1280 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1282 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1283 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1286 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1287 write (iout,'(a,f10.4,4(16x,f10.4))') &
1288 'Center ',(bsc(j,i),j=1,nlobi)
1289 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1298 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1299 ! added by Urszula Kozlowska 07/11/2007
1301 !el Maximum number of SC local term fitting function coefficiants
1302 !el integer,parameter :: maxsccoef=65
1304 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1307 read (irotam,*,end=112,err=112)
1309 read (irotam,*,end=112,err=112)
1312 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1317 ! Read the parameters of the probability distribution/energy expression
1318 ! of the side chains.
1320 write (2,*) "Start reading ROTAM_PDB"
1322 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1326 dsc_inv(i)=1.0D0/dsc(i)
1337 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1338 ((blower(k,l,1),l=1,k),k=1,3)
1340 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1341 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1342 ((blower(k,l,j),l=1,k),k=1,3)
1349 akl=akl+blower(k,m,j)*blower(l,m,j)
1359 write (2,*) "End reading ROTAM_PDB"
1365 ! Read torsional parameters in old format
1367 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1369 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1370 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1371 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1373 !el from energy module--------
1374 allocate(v1(nterm_old,ntortyp,ntortyp))
1375 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1376 !el---------------------------
1381 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1387 write (iout,'(/a/)') 'Torsional constants:'
1390 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1391 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1397 ! Read torsional parameters
1399 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1400 read (itorp,*,end=113,err=113) ntortyp
1401 !el from energy module---------
1402 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1403 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1405 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1406 allocate(vlor2(maxlor,ntortyp,ntortyp))
1407 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1408 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1410 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1411 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1412 !el---------------------------
1415 !el---------------------------
1417 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1419 itortyp(i)=-itortyp(-i)
1421 ! itortyp(ntyp1)=ntortyp
1422 ! itortyp(-ntyp1)=-ntortyp
1424 write (iout,*) 'ntortyp',ntortyp
1426 do j=-ntortyp+1,ntortyp-1
1427 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1429 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1430 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1433 do k=1,nterm(i,j,iblock)
1434 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1436 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1437 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1438 v0ij=v0ij+si*v1(k,i,j,iblock)
1440 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1441 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1442 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1444 do k=1,nlor(i,j,iblock)
1445 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1446 vlor2(k,i,j),vlor3(k,i,j)
1447 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1450 v0(-i,-j,iblock)=v0ij
1456 write (iout,'(/a/)') 'Torsional constants:'
1458 do i=-ntortyp,ntortyp
1459 do j=-ntortyp,ntortyp
1460 write (iout,*) 'ityp',i,' jtyp',j
1461 write (iout,*) 'Fourier constants'
1462 do k=1,nterm(i,j,iblock)
1463 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1466 write (iout,*) 'Lorenz constants'
1467 do k=1,nlor(i,j,iblock)
1468 write (iout,'(3(1pe15.5))') &
1469 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1475 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1477 ! 6/23/01 Read parameters for double torsionals
1479 !el from energy module------------
1480 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1481 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1482 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1483 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1484 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1485 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1486 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1487 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1488 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1489 !---------------------------------
1493 do j=-ntortyp+1,ntortyp-1
1494 do k=-ntortyp+1,ntortyp-1
1495 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1496 ! write (iout,*) "OK onelett",
1499 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1500 .or. t3.ne.toronelet(k)) then
1501 write (iout,*) "Error in double torsional parameter file",&
1504 call MPI_Finalize(Ierror)
1506 stop "Error in double torsional parameter file"
1508 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1509 ntermd_2(i,j,k,iblock)
1510 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1511 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1512 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1513 ntermd_1(i,j,k,iblock))
1514 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1515 ntermd_1(i,j,k,iblock))
1516 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1517 ntermd_1(i,j,k,iblock))
1518 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1519 ntermd_1(i,j,k,iblock))
1520 ! Martix of D parameters for one dimesional foureir series
1521 do l=1,ntermd_1(i,j,k,iblock)
1522 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1523 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1524 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1525 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1526 ! write(iout,*) "whcodze" ,
1527 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1529 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1530 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1531 v2s(m,l,i,j,k,iblock),&
1532 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1533 ! Martix of D parameters for two dimesional fourier series
1534 do l=1,ntermd_2(i,j,k,iblock)
1536 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1537 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1538 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1539 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1548 write (iout,*) 'Constants for double torsionals'
1551 do j=-ntortyp+1,ntortyp-1
1552 do k=-ntortyp+1,ntortyp-1
1553 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1554 ' nsingle',ntermd_1(i,j,k,iblock),&
1555 ' ndouble',ntermd_2(i,j,k,iblock)
1557 write (iout,*) 'Single angles:'
1558 do l=1,ntermd_1(i,j,k,iblock)
1559 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1560 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1561 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1562 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1565 write (iout,*) 'Pairs of angles:'
1566 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1567 do l=1,ntermd_2(i,j,k,iblock)
1568 write (iout,'(i5,20f10.5)') &
1569 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1572 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1573 do l=1,ntermd_2(i,j,k,iblock)
1574 write (iout,'(i5,20f10.5)') &
1575 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1576 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1585 ! Read of Side-chain backbone correlation parameters
1586 ! Modified 11 May 2012 by Adasko
1589 read (isccor,*,end=119,err=119) nsccortyp
1591 !el from module energy-------------
1592 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1593 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1594 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1595 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1596 !-----------------------------------
1598 !el from module energy-------------
1599 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1601 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1603 isccortyp(i)=-isccortyp(-i)
1605 iscprol=isccortyp(20)
1606 ! write (iout,*) 'ntortyp',ntortyp
1608 !c maxinter is maximum interaction sites
1609 !el from module energy---------
1610 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1611 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1612 -nsccortyp:nsccortyp))
1613 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1614 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1615 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1616 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1617 !-----------------------------------
1621 read (isccor,*,end=119,err=119) &
1622 nterm_sccor(i,j),nlor_sccor(i,j)
1628 nterm_sccor(-i,j)=nterm_sccor(i,j)
1629 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1630 nterm_sccor(i,-j)=nterm_sccor(i,j)
1631 do k=1,nterm_sccor(i,j)
1632 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1634 if (j.eq.iscprol) then
1635 if (i.eq.isccortyp(10)) then
1636 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1637 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1639 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1640 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1641 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1642 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1643 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1644 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1645 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1646 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1649 if (i.eq.isccortyp(10)) then
1650 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1651 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1653 if (j.eq.isccortyp(10)) then
1654 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1655 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1657 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1658 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1659 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1660 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1661 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1662 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1666 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1667 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1668 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1669 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1672 do k=1,nlor_sccor(i,j)
1673 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1674 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1675 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1676 (1+vlor3sccor(k,i,j)**2)
1678 v0sccor(l,i,j)=v0ijsccor
1679 v0sccor(l,-i,j)=v0ijsccor1
1680 v0sccor(l,i,-j)=v0ijsccor2
1681 v0sccor(l,-i,-j)=v0ijsccor3
1687 !el from module energy-------------
1688 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1690 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1691 ! write (iout,*) 'ntortyp',ntortyp
1693 !c maxinter is maximum interaction sites
1694 !el from module energy---------
1695 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1696 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1697 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1698 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1699 !-----------------------------------
1703 read (isccor,*,end=119,err=119) &
1704 nterm_sccor(i,j),nlor_sccor(i,j)
1708 do k=1,nterm_sccor(i,j)
1709 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1711 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1714 do k=1,nlor_sccor(i,j)
1715 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1716 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1717 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1718 (1+vlor3sccor(k,i,j)**2)
1720 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1728 write (iout,'(/a/)') 'Torsional constants:'
1731 write (iout,*) 'ityp',i,' jtyp',j
1732 write (iout,*) 'Fourier constants'
1733 do k=1,nterm_sccor(i,j)
1734 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1736 write (iout,*) 'Lorenz constants'
1737 do k=1,nlor_sccor(i,j)
1738 write (iout,'(3(1pe15.5))') &
1739 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1746 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1747 ! interaction energy of the Gly, Ala, and Pro prototypes.
1752 write (iout,*) "Coefficients of the cumulants"
1754 read (ifourier,*) nloctyp
1755 !write(iout,*) "nloctyp",nloctyp
1756 !el from module energy-------
1757 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1758 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1759 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1760 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1761 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1762 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1763 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1764 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1768 !--------------------------------
1771 read (ifourier,*,end=115,err=115)
1772 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1774 write (iout,*) 'Type',i
1775 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1785 B1tilde(1,-i) =-b(3)
1787 ! b1tilde(1,i)=0.0d0
1788 ! b1tilde(2,i)=0.0d0
1813 Ctilde(1,2,-i)=-b(9)
1817 ! Ctilde(1,1,i)=0.0d0
1818 ! Ctilde(1,2,i)=0.0d0
1819 ! Ctilde(2,1,i)=0.0d0
1820 ! Ctilde(2,2,i)=0.0d0
1838 Dtilde(1,2,-i)=-b(8)
1842 ! Dtilde(1,1,i)=0.0d0
1843 ! Dtilde(1,2,i)=0.0d0
1844 ! Dtilde(2,1,i)=0.0d0
1845 ! Dtilde(2,2,i)=0.0d0
1846 EE(1,1,i)= b(10)+b(11)
1847 EE(2,2,i)=-b(10)+b(11)
1848 EE(2,1,i)= b(12)-b(13)
1849 EE(1,2,i)= b(12)+b(13)
1850 EE(1,1,-i)= b(10)+b(11)
1851 EE(2,2,-i)=-b(10)+b(11)
1852 EE(2,1,-i)=-b(12)+b(13)
1853 EE(1,2,-i)=-b(12)-b(13)
1859 ! ee(2,1,i)=ee(1,2,i)
1863 write (iout,*) 'Type',i
1865 write(iout,*) B1(1,i),B1(2,i)
1867 write(iout,*) B2(1,i),B2(2,i)
1870 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1874 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1878 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1883 ! Read electrostatic-interaction parameters
1888 write (iout,'(/a)') 'Electrostatic interaction constants:'
1889 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1890 'IT','JT','APP','BPP','AEL6','AEL3'
1892 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1893 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1894 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1895 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1900 app (i,j)=epp(i,j)*rri*rri
1901 bpp (i,j)=-2.0D0*epp(i,j)*rri
1902 ael6(i,j)=elpp6(i,j)*4.2D0**6
1903 ael3(i,j)=elpp3(i,j)*4.2D0**3
1905 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1911 ! Read side-chain interaction parameters.
1913 !el from module energy - COMMON.INTERACT-------
1914 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1915 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1916 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1917 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1918 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1927 !--------------------------------
1929 read (isidep,*,end=117,err=117) ipot,expon
1930 if (ipot.lt.1 .or. ipot.gt.5) then
1931 write (iout,'(2a)') 'Error while reading SC interaction',&
1932 'potential file - unknown potential type.'
1934 call MPI_Finalize(Ierror)
1940 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1941 ', exponents are ',expon,2*expon
1942 ! goto (10,20,30,30,40) ipot
1944 !----------------------- LJ potential ---------------------------------
1946 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1947 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1948 (sigma0(i),i=1,ntyp)
1950 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1951 write (iout,'(a/)') 'The epsilon array:'
1952 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1953 write (iout,'(/a)') 'One-body parameters:'
1954 write (iout,'(a,4x,a)') 'residue','sigma'
1955 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1958 !----------------------- LJK potential --------------------------------
1960 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1962 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1964 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1965 write (iout,'(a/)') 'The epsilon array:'
1966 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1967 write (iout,'(/a)') 'One-body parameters:'
1968 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1969 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1973 !---------------------- GB or BP potential -----------------------------
1977 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1979 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1980 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1981 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1982 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1983 ! For the GB potential convert sigma'**2 into chi'
1986 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1990 write (iout,'(/a/)') 'Parameters of the BP potential:'
1991 write (iout,'(a/)') 'The epsilon array:'
1992 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1993 write (iout,'(/a)') 'One-body parameters:'
1994 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1996 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1997 chip(i),alp(i),i=1,ntyp)
2000 !--------------------- GBV potential -----------------------------------
2002 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2003 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2004 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2005 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2007 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2008 write (iout,'(a/)') 'The epsilon array:'
2009 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010 write (iout,'(/a)') 'One-body parameters:'
2011 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2012 's||/s_|_^2',' chip ',' alph '
2013 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2014 sigii(i),chip(i),alp(i),i=1,ntyp)
2017 write(iout,*)"Wrong ipot"
2022 !-----------------------------------------------------------------------
2023 ! Calculate the "working" parameters of SC interactions.
2025 !el from module energy - COMMON.INTERACT-------
2026 allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2027 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2034 !--------------------------------
2043 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2044 sigma(j,i)=sigma(i,j)
2045 rs0(i,j)=dwa16*sigma(i,j)
2049 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2050 'Working parameters of the SC interactions:',&
2051 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2056 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2065 sigeps=dsign(1.0D0,epsij)
2067 aa(i,j)=epsij*rrij*rrij
2068 bb(i,j)=-sigeps*epsij*rrij
2072 sigt1sq=sigma0(i)**2
2073 sigt2sq=sigma0(j)**2
2076 ratsig1=sigt2sq/sigt1sq
2077 ratsig2=1.0D0/ratsig1
2078 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2079 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2080 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2084 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2085 sigmaii(i,j)=rsum_max
2086 sigmaii(j,i)=rsum_max
2088 ! sigmaii(i,j)=r0(i,j)
2089 ! sigmaii(j,i)=r0(i,j)
2091 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2092 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2093 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2094 augm(i,j)=epsij*r_augm**(2*expon)
2095 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2102 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2103 restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
2104 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2109 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2114 ! Define the SC-p interaction constants (hard-coded; old style)
2117 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2119 ! aad(i,1)=0.3D0*4.0D0**12
2120 ! Following line for constants currently implemented
2121 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2122 aad(i,1)=1.5D0*4.0D0**12
2123 ! aad(i,1)=0.17D0*5.6D0**12
2125 ! "Soft" SC-p repulsion
2127 ! Following line for constants currently implemented
2128 ! aad(i,1)=0.3D0*4.0D0**6
2129 ! "Hard" SC-p repulsion
2130 bad(i,1)=3.0D0*4.0D0**6
2131 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2140 ! 8/9/01 Read the SC-p interaction constants from file
2143 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2146 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2147 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2148 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2149 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2153 write (iout,*) "Parameters of SC-p interactions:"
2155 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2156 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2162 ! Define the constants of the disulfide bridge
2166 ! Old arbitrary potential - commented out.
2171 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2172 ! energy surface of diethyl disulfide.
2173 ! A. Liwo and U. Kozlowska, 11/24/03
2190 write (iout,'(/a)') "Disulfide bridge parameters:"
2191 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2192 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2193 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2194 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2198 111 write (iout,*) "Error reading bending energy parameters."
2200 112 write (iout,*) "Error reading rotamer energy parameters."
2202 113 write (iout,*) "Error reading torsional energy parameters."
2204 114 write (iout,*) "Error reading double torsional energy parameters."
2206 115 write (iout,*) &
2207 "Error reading cumulant (multibody energy) parameters."
2209 116 write (iout,*) "Error reading electrostatic energy parameters."
2211 117 write (iout,*) "Error reading side chain interaction parameters."
2213 118 write (iout,*) "Error reading SCp interaction parameters."
2215 119 write (iout,*) "Error reading SCCOR parameters"
2218 call MPI_Finalize(Ierror)
2222 end subroutine parmread
2224 !-----------------------------------------------------------------------------
2226 !-----------------------------------------------------------------------------
2227 subroutine printmat(ldim,m,n,iout,key,a)
2230 character(len=3),dimension(n) :: key
2231 real(kind=8),dimension(ldim,n) :: a
2233 integer :: i,j,k,m,iout,nlim
2237 write (iout,1000) (key(k),k=i,nlim)
2239 1000 format (/5x,8(6x,a3))
2240 1020 format (/80(1h-)/)
2242 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2245 1010 format (a3,2x,8(f9.4))
2247 end subroutine printmat
2248 !-----------------------------------------------------------------------------
2250 !-----------------------------------------------------------------------------
2252 ! Read the PDB file and convert the peptide geometry into virtual-chain
2255 use energy_data, only: itype
2259 use control, only: rescode
2260 ! implicit real*8 (a-h,o-z)
2261 ! include 'DIMENSIONS'
2262 ! include 'COMMON.LOCAL'
2263 ! include 'COMMON.VAR'
2264 ! include 'COMMON.CHAIN'
2265 ! include 'COMMON.INTERACT'
2266 ! include 'COMMON.IOUNITS'
2267 ! include 'COMMON.GEO'
2268 ! include 'COMMON.NAMES'
2269 ! include 'COMMON.CONTROL'
2270 ! include 'COMMON.DISTFIT'
2271 ! include 'COMMON.SETUP'
2272 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2274 logical :: lprn=.true.,fail
2275 real(kind=8),dimension(3) :: e1,e2,e3
2276 real(kind=8) :: dcj,efree_temp
2277 character(len=3) :: seq,res
2278 character(len=5) :: atom
2279 character(len=80) :: card
2280 real(kind=8),dimension(3,20) :: sccor
2281 integer :: kkk,lll,icha,kupa !rescode,
2284 integer,dimension(2,maxres/3) :: hfrag_alloc
2285 integer,dimension(4,maxres/3) :: bfrag_alloc
2286 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2292 ! write (2,*) "UNRES_PDB",unres_pdb
2300 !-----------------------------
2301 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2302 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2305 read (ipdbin,'(a80)',end=10) card
2306 ! write (iout,'(a)') card
2307 if (card(:5).eq.'HELIX') then
2310 read(card(22:25),*) hfrag(1,nhfrag)
2311 read(card(34:37),*) hfrag(2,nhfrag)
2313 if (card(:5).eq.'SHEET') then
2316 read(card(24:26),*) bfrag(1,nbfrag)
2317 read(card(35:37),*) bfrag(2,nbfrag)
2318 !rc----------------------------------------
2319 !rc to be corrected !!!
2320 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2321 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2322 !rc----------------------------------------
2324 if (card(:3).eq.'END') then
2326 else if (card(:3).eq.'TER') then
2330 itype(ires_old)=ntyp1
2332 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2335 dc(j,ires)=sccor(j,iii)
2338 call sccenter(ires,iii,sccor)
2343 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2344 ! Fish out the ATOM cards.
2345 if (index(card(1:4),'ATOM').gt.0) then
2346 read (card(12:16),*) atom
2347 ! write (iout,*) "! ",atom," !",ires
2348 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2349 read (card(23:26),*) ires
2350 read (card(18:20),'(a3)') res
2351 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2352 ! & " ires_old",ires_old
2353 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2354 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2355 if (ires-ishift+ishift1.ne.ires_old) then
2356 ! Calculate the CM of the preceding residue.
2357 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2359 ! write (iout,*) "Calculating sidechain center iii",iii
2362 dc(j,ires+nres)=sccor(j,iii)
2365 call sccenter(ires_old,iii,sccor)
2369 ! Start new residue.
2370 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2373 else if (ibeg.eq.1) then
2374 write (iout,*) "BEG ires",ires
2376 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2380 ires=ires-ishift+ishift1
2382 ! write (iout,*) "ishift",ishift," ires",ires,&
2383 ! " ires_old",ires_old
2385 else if (ibeg.eq.2) then
2387 ishift=-ires_old+ires-1 !!!!!
2388 ishift1=ishift1-1 !!!!!
2389 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2390 ires=ires-ishift+ishift1
2394 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2395 ires=ires-ishift+ishift1
2398 if (res.eq.'ACE' .or. res.eq.'NHE') then
2401 itype(ires)=rescode(ires,res,0)
2404 ires=ires-ishift+ishift1
2406 ! write (iout,*) "ires_old",ires_old," ires",ires
2407 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2410 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2411 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2412 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2413 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2414 ! write (iout,*) "backbone ",atom
2416 write (iout,'(2i3,2x,a,3f8.3)') &
2417 ires,itype(ires),res,(c(j,ires),j=1,3)
2421 sccor(j,iii)=c(j,ires)
2423 ! write (*,*) card(23:27),ires,itype(ires)
2424 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2425 atom.ne.'N' .and. atom.ne.'C' .and. &
2426 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2427 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2428 ! write (iout,*) "sidechain ",atom
2430 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2434 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2435 if (ires.eq.0) return
2436 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2440 ! write (iout,*) i,itype(i)
2441 if (itype(i).eq.ntyp1) then
2442 ! write (iout,*) "dummy",i,itype(i)
2444 c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2445 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2450 ! Calculate the CM of the last side chain.
2454 dc(j,ires)=sccor(j,iii)
2457 call sccenter(ires,iii,sccor)
2463 if (itype(nres).ne.10) then
2467 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2468 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2475 c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
2479 dcj=c(j,nres-2)-c(j,nres-3)
2480 c(j,nres)=c(j,nres-1)+dcj
2481 c(j,2*nres)=c(j,nres)
2485 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2487 if (nres.ne.nres0) then
2488 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2490 stop "Error nres value in WHAM input"
2493 !---------------------------------
2494 !el reallocate tables
2497 ! hfrag_alloc(j,i)=hfrag(j,i)
2500 ! bfrag_alloc(j,i)=bfrag(j,i)
2506 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2507 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2508 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2509 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2513 ! hfrag(j,i)=hfrag_alloc(j,i)
2518 ! bfrag(j,i)=bfrag_alloc(j,i)
2521 !el end reallocate tables
2522 !---------------------------------
2530 c(j,2*nres)=c(j,nres)
2532 if (itype(1).eq.ntyp1) then
2536 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2537 call refsys(2,3,4,e1,e2,e3,fail)
2544 c(j,1)=c(j,2)-3.8d0*e2(j)
2554 ! Copy the coordinates to reference coordinates
2560 ! Calculate internal coordinates.
2562 write (iout,'(/a)') &
2563 "Cartesian coordinates of the reference structure"
2564 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2565 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2567 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2568 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2569 (c(j,ires+nres),j=1,3)
2572 ! znamy już nres wiec mozna alokowac tablice
2573 ! Calculate internal coordinates.
2574 if(me.eq.king.or..not.out1file)then
2575 write (iout,'(a)') &
2576 "Backbone and SC coordinates as read from the PDB"
2578 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2579 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2580 (c(j,nres+ires),j=1,3)
2584 if(.not.allocated(vbld)) then
2585 allocate(vbld(2*nres))
2590 if(.not.allocated(vbld_inv)) then
2591 allocate(vbld_inv(2*nres))
2597 if(.not.allocated(theta)) then
2598 allocate(theta(nres+2))
2602 if(.not.allocated(phi)) allocate(phi(nres+2))
2603 if(.not.allocated(alph)) allocate(alph(nres+2))
2604 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2605 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2606 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2607 if(.not.allocated(costtab)) allocate(costtab(nres))
2608 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2609 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2610 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2611 if(.not.allocated(xxref)) allocate(xxref(nres))
2612 if(.not.allocated(yyref)) allocate(yyref(nres))
2613 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2614 if(.not.allocated(dc_norm)) then
2615 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2616 allocate(dc_norm(3,0:2*nres+2))
2620 call int_from_cart(.true.,.false.)
2621 call sc_loc_geom(.false.)
2623 thetaref(i)=theta(i)
2633 dc(j,i)=c(j,i+1)-c(j,i)
2634 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2639 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2640 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2642 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2646 ! Copy the coordinates to reference coordinates
2647 ! Splits to single chain if occurs
2651 ! cref(j,i,cou)=c(j,i)
2655 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2656 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2657 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2658 !-----------------------------
2664 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2666 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2669 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2674 cref(j,i,cou)=c(j,i)
2675 cref(j,i+nres,cou)=c(j,i+nres)
2677 chain_rep(j,lll,kkk)=c(j,i)
2678 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2682 write (iout,*) chain_length
2683 if (chain_length.eq.0) chain_length=nres
2685 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2686 chain_rep(j,chain_length+nres,symetr) &
2687 =chain_rep(j,chain_length+nres,1)
2690 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2692 ! do kkk=1,chain_length
2693 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2697 ! makes copy of chains
2698 write (iout,*) "symetr", symetr
2700 if (symetr.gt.1) then
2707 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2713 ! write (iout,*) i,icha
2714 do lll=1,chain_length
2716 if (cou.le.nres) then
2718 kupa=mod(lll,chain_length)
2719 iprzes=(kkk-1)*chain_length+lll
2720 if (kupa.eq.0) kupa=chain_length
2721 ! write (iout,*) "kupa", kupa
2722 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2723 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2730 !-koniec robienia kopii
2733 write (iout,*) "nowa struktura", nperm
2735 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2737 cref(3,i,kkk),cref(1,nres+i,kkk),&
2738 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2740 100 format (//' alpha-carbon coordinates ',&
2741 ' centroid coordinates'/ &
2742 ' ', 6X,'X',11X,'Y',11X,'Z', &
2743 10X,'X',11X,'Y',11X,'Z')
2744 110 format (a,'(',i3,')',6f12.5)
2750 bfrag(i,j)=bfrag(i,j)-ishift
2756 hfrag(i,j)=hfrag(i,j)-ishift
2762 end subroutine readpdb
2763 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2764 !-----------------------------------------------------------------------------
2766 !-----------------------------------------------------------------------------
2767 subroutine read_control
2781 use random, only: random_init
2782 ! implicit real*8 (a-h,o-z)
2783 ! include 'DIMENSIONS'
2785 use prng, only:prng_restart
2787 logical :: OKRandom!, prng_restart
2790 ! include 'COMMON.IOUNITS'
2791 ! include 'COMMON.TIME1'
2792 ! include 'COMMON.THREAD'
2793 ! include 'COMMON.SBRIDGE'
2794 ! include 'COMMON.CONTROL'
2795 ! include 'COMMON.MCM'
2796 ! include 'COMMON.MAP'
2797 ! include 'COMMON.HEADER'
2798 ! include 'COMMON.CSA'
2799 ! include 'COMMON.CHAIN'
2800 ! include 'COMMON.MUCA'
2801 ! include 'COMMON.MD'
2802 ! include 'COMMON.FFIELD'
2803 ! include 'COMMON.INTERACT'
2804 ! include 'COMMON.SETUP'
2805 !el integer :: KDIAG,ICORFL,IXDR
2806 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2807 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2808 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2809 ! character(len=80) :: ucase
2810 character(len=640) :: controlcard
2812 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2818 read (INP,'(a)') titel
2819 call card_concat(controlcard,.true.)
2820 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2821 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2822 call reada(controlcard,'SEED',seed,0.0D0)
2823 call random_init(seed)
2824 ! Set up the time limit (caution! The time must be input in minutes!)
2825 read_cart=index(controlcard,'READ_CART').gt.0
2826 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2827 call readi(controlcard,'SYM',symetr,1)
2828 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2829 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2830 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2831 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2832 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2833 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2834 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2835 call reada(controlcard,'DRMS',drms,0.1D0)
2836 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2837 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2838 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2839 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2840 write (iout,'(a,f10.1)')'DRMS = ',drms
2841 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2842 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2844 call readi(controlcard,'NZ_START',nz_start,0)
2845 call readi(controlcard,'NZ_END',nz_end,0)
2846 call readi(controlcard,'IZ_SC',iz_sc,0)
2847 timlim=60.0D0*timlim
2848 safety = 60.0d0*safety
2851 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2852 minim=(index(controlcard,'MINIMIZE').gt.0)
2853 dccart=(index(controlcard,'CART').gt.0)
2854 overlapsc=(index(controlcard,'OVERLAP').gt.0)
2855 overlapsc=.not.overlapsc
2856 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2857 searchsc=.not.searchsc
2858 sideadd=(index(controlcard,'SIDEADD').gt.0)
2859 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2860 outpdb=(index(controlcard,'PDBOUT').gt.0)
2861 outmol2=(index(controlcard,'MOL2OUT').gt.0)
2862 pdbref=(index(controlcard,'PDBREF').gt.0)
2863 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2864 indpdb=index(controlcard,'PDBSTART')
2865 extconf=(index(controlcard,'EXTCONF').gt.0)
2866 call readi(controlcard,'IPRINT',iprint,0)
2867 call readi(controlcard,'MAXGEN',maxgen,10000)
2868 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2869 call readi(controlcard,"KDIAG",kdiag,0)
2870 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2871 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2872 write (iout,*) "RESCALE_MODE",rescale_mode
2873 split_ene=index(controlcard,'SPLIT_ENE').gt.0
2874 if (index(controlcard,'REGULAR').gt.0.0D0) then
2875 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2879 if (index(controlcard,'CHECKGRAD').gt.0) then
2881 if (index(controlcard,'CART').gt.0) then
2883 elseif (index(controlcard,'CARINT').gt.0) then
2888 elseif (index(controlcard,'THREAD').gt.0) then
2890 call readi(controlcard,'THREAD',nthread,0)
2891 if (nthread.gt.0) then
2892 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2895 write (iout,'(a)')'A number has to follow the THREAD keyword.'
2896 stop 'Error termination in Read_Control.'
2898 else if (index(controlcard,'MCMA').gt.0) then
2900 else if (index(controlcard,'MCEE').gt.0) then
2902 else if (index(controlcard,'MULTCONF').gt.0) then
2904 else if (index(controlcard,'MAP').gt.0) then
2906 call readi(controlcard,'MAP',nmap,0)
2907 else if (index(controlcard,'CSA').gt.0) then
2909 !rc else if (index(controlcard,'ZSCORE').gt.0) then
2911 !rc ZSCORE is rm from UNRES, modecalc=9 is available
2914 !fcm else if (index(controlcard,'MCMF').gt.0) then
2916 else if (index(controlcard,'SOFTREG').gt.0) then
2918 else if (index(controlcard,'CHECK_BOND').gt.0) then
2920 else if (index(controlcard,'TEST').gt.0) then
2922 else if (index(controlcard,'MD').gt.0) then
2924 else if (index(controlcard,'RE ').gt.0) then
2928 lmuca=index(controlcard,'MUCA').gt.0
2929 call readi(controlcard,'MUCADYN',mucadyn,0)
2930 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
2931 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
2933 write (iout,*) 'MUCADYN=',mucadyn
2934 write (iout,*) 'MUCASMOOTH=',muca_smooth
2937 iscode=index(controlcard,'ONE_LETTER')
2938 indphi=index(controlcard,'PHI')
2939 indback=index(controlcard,'BACK')
2940 iranconf=index(controlcard,'RAND_CONF')
2941 i2ndstr=index(controlcard,'USE_SEC_PRED')
2942 gradout=index(controlcard,'GRADOUT').gt.0
2943 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
2944 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
2945 if (me.eq.king .or. .not.out1file ) &
2946 write (iout,*) "DISTCHAINMAX",distchainmax
2948 if(me.eq.king.or..not.out1file) &
2949 write (iout,'(2a)') diagmeth(kdiag),&
2950 ' routine used to diagonalize matrices.'
2952 end subroutine read_control
2953 !-----------------------------------------------------------------------------
2954 subroutine read_REMDpar
2956 ! Read REMD settings
2963 use control_data, only:out1file
2964 ! implicit real*8 (a-h,o-z)
2965 ! include 'DIMENSIONS'
2966 ! include 'COMMON.IOUNITS'
2967 ! include 'COMMON.TIME1'
2968 ! include 'COMMON.MD'
2971 !el include 'COMMON.LANGEVIN'
2973 !el include 'COMMON.LANGEVIN.lang0'
2975 ! include 'COMMON.INTERACT'
2976 ! include 'COMMON.NAMES'
2977 ! include 'COMMON.GEO'
2978 ! include 'COMMON.REMD'
2979 ! include 'COMMON.CONTROL'
2980 ! include 'COMMON.SETUP'
2981 ! character(len=80) :: ucase
2982 character(len=320) :: controlcard
2983 character(len=3200) :: controlcard1
2984 integer :: iremd_m_total
2987 ! real(kind=8) :: var,ene
2989 if(me.eq.king.or..not.out1file) &
2990 write (iout,*) "REMD setup"
2992 call card_concat(controlcard,.true.)
2993 call readi(controlcard,"NREP",nrep,3)
2994 call readi(controlcard,"NSTEX",nstex,1000)
2995 call reada(controlcard,"RETMIN",retmin,10.0d0)
2996 call reada(controlcard,"RETMAX",retmax,1000.0d0)
2997 mremdsync=(index(controlcard,'SYNC').gt.0)
2998 call readi(controlcard,"NSYN",i_sync_step,100)
2999 restart1file=(index(controlcard,'REST1FILE').gt.0)
3000 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3001 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3002 if(max_cache_traj_use.gt.max_cache_traj) &
3003 max_cache_traj_use=max_cache_traj
3004 if(me.eq.king.or..not.out1file) then
3005 !d if (traj1file) then
3006 !rc caching is in testing - NTWX is not ignored
3007 !d write (iout,*) "NTWX value is ignored"
3008 !d write (iout,*) " trajectory is stored to one file by master"
3009 !d write (iout,*) " before exchange at NSTEX intervals"
3011 write (iout,*) "NREP= ",nrep
3012 write (iout,*) "NSTEX= ",nstex
3013 write (iout,*) "SYNC= ",mremdsync
3014 write (iout,*) "NSYN= ",i_sync_step
3015 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3018 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3019 if (index(controlcard,'TLIST').gt.0) then
3021 call card_concat(controlcard1,.true.)
3022 read(controlcard1,*) (remd_t(i),i=1,nrep)
3023 if(me.eq.king.or..not.out1file) &
3024 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3027 if (index(controlcard,'MLIST').gt.0) then
3029 call card_concat(controlcard1,.true.)
3030 read(controlcard1,*) (remd_m(i),i=1,nrep)
3031 if(me.eq.king.or..not.out1file) then
3032 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3035 iremd_m_total=iremd_m_total+remd_m(i)
3037 write (iout,*) 'Total number of replicas ',iremd_m_total
3040 if(me.eq.king.or..not.out1file) &
3041 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3043 end subroutine read_REMDpar
3044 !-----------------------------------------------------------------------------
3045 subroutine read_MDpar
3049 use control_data, only: r_cut,rlamb,out1file
3051 use geometry_data, only: pi
3053 ! implicit real*8 (a-h,o-z)
3054 ! include 'DIMENSIONS'
3055 ! include 'COMMON.IOUNITS'
3056 ! include 'COMMON.TIME1'
3057 ! include 'COMMON.MD'
3060 !el include 'COMMON.LANGEVIN'
3062 !el include 'COMMON.LANGEVIN.lang0'
3064 ! include 'COMMON.INTERACT'
3065 ! include 'COMMON.NAMES'
3066 ! include 'COMMON.GEO'
3067 ! include 'COMMON.SETUP'
3068 ! include 'COMMON.CONTROL'
3069 ! include 'COMMON.SPLITELE'
3070 ! character(len=80) :: ucase
3071 character(len=320) :: controlcard
3076 call card_concat(controlcard,.true.)
3077 call readi(controlcard,"NSTEP",n_timestep,1000000)
3078 call readi(controlcard,"NTWE",ntwe,100)
3079 call readi(controlcard,"NTWX",ntwx,1000)
3080 call reada(controlcard,"DT",d_time,1.0d-1)
3081 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3082 call reada(controlcard,"DAMAX",damax,1.0d1)
3083 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3084 call readi(controlcard,"LANG",lang,0)
3085 RESPA = index(controlcard,"RESPA") .gt. 0
3086 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3087 ntime_split0=ntime_split
3088 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3089 ntime_split0=ntime_split
3090 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3091 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3092 rest = index(controlcard,"REST").gt.0
3093 tbf = index(controlcard,"TBF").gt.0
3094 usampl = index(controlcard,"USAMPL").gt.0
3095 mdpdb = index(controlcard,"MDPDB").gt.0
3096 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3097 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3098 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3099 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3100 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3101 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3102 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3103 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3104 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3105 large = index(controlcard,"LARGE").gt.0
3106 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3107 rattle = index(controlcard,"RATTLE").gt.0
3108 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3114 if(me.eq.king.or..not.out1file) then
3116 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3118 write (iout,'(a)') "The units are:"
3119 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3120 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3121 " acceleration: angstrom/(48.9 fs)**2"
3122 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3124 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3125 write (iout,'(a60,f10.5,a)') &
3126 "Initial time step of numerical integration:",d_time,&
3128 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3130 write (iout,'(2a,i4,a)') &
3131 "A-MTS algorithm used; initial time step for fast-varying",&
3132 " short-range forces split into",ntime_split," steps."
3133 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3134 r_cut," lambda",rlamb
3136 write (iout,'(2a,f10.5)') &
3137 "Maximum acceleration threshold to reduce the time step",&
3138 "/increase split number:",damax
3139 write (iout,'(2a,f10.5)') &
3140 "Maximum predicted energy drift to reduce the timestep",&
3141 "/increase split number:",edriftmax
3142 write (iout,'(a60,f10.5)') &
3143 "Maximum velocity threshold to reduce velocities:",dvmax
3144 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3145 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3146 if (rattle) write (iout,'(a60)') &
3147 "Rattle algorithm used to constrain the virtual bonds"
3151 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3152 call reada(controlcard,"RWAT",rwat,1.4d0)
3153 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3154 surfarea=index(controlcard,"SURFAREA").gt.0
3155 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3156 if(me.eq.king.or..not.out1file)then
3157 write (iout,'(/a,$)') "Langevin dynamics calculation"
3159 write (iout,'(a/)') &
3160 " with direct integration of Langevin equations"
3161 else if (lang.eq.2) then
3162 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3163 else if (lang.eq.3) then
3164 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3165 else if (lang.eq.4) then
3166 write (iout,'(a/)') " in overdamped mode"
3168 write (iout,'(//a,i5)') &
3169 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3172 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3173 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3174 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3175 write (iout,'(a60,f10.5)') &
3176 "Scaling factor of the friction forces:",scal_fric
3177 if (surfarea) write (iout,'(2a,i10,a)') &
3178 "Friction coefficients will be scaled by solvent-accessible",&
3179 " surface area every",reset_fricmat," steps."
3181 ! Calculate friction coefficients and bounds of stochastic forces
3182 eta=6*pi*cPoise*etawat
3183 if(me.eq.king.or..not.out1file) &
3184 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3186 gamp=scal_fric*(pstok+rwat)*eta
3187 stdfp=dsqrt(2*Rb*t_bath/d_time)
3188 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3190 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3191 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3193 if(me.eq.king.or..not.out1file)then
3194 write (iout,'(/2a/)') &
3195 "Radii of site types and friction coefficients and std's of",&
3196 " stochastic forces of fully exposed sites"
3197 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3199 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3200 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3204 if(me.eq.king.or..not.out1file)then
3205 write (iout,'(a)') "Berendsen bath calculation"
3206 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3207 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3209 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3210 count_reset_moment," steps"
3212 write (iout,'(a,i10,a)') &
3213 "Velocities will be reset at random every",count_reset_vel,&
3217 if(me.eq.king.or..not.out1file) &
3218 write (iout,'(a31)') "Microcanonical mode calculation"
3220 if(me.eq.king.or..not.out1file)then
3221 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3223 write(iout,*) "MD running with constraints."
3224 write(iout,*) "Equilibration time ", eq_time, " mtus."
3225 write(iout,*) "Constraining ", nfrag," fragments."
3226 write(iout,*) "Length of each fragment, weight and q0:"
3228 write (iout,*) "Set of restraints #",iset
3230 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3231 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3233 write(iout,*) "constraints between ", npair, "fragments."
3234 write(iout,*) "constraint pairs, weights and q0:"
3236 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3237 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3239 write(iout,*) "angle constraints within ", nfrag_back,&
3240 "backbone fragments."
3241 write(iout,*) "fragment, weights:"
3243 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3244 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3245 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3248 iset=mod(kolor,nset)+1
3251 if(me.eq.king.or..not.out1file) &
3252 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3254 end subroutine read_MDpar
3255 !-----------------------------------------------------------------------------
3259 ! implicit real*8 (a-h,o-z)
3260 ! include 'DIMENSIONS'
3261 ! include 'COMMON.MAP'
3262 ! include 'COMMON.IOUNITS'
3263 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3264 character(len=80) :: mapcard !,ucase
3267 ! real(kind=8) :: var,ene
3270 read (inp,'(a)') mapcard
3271 mapcard=ucase(mapcard)
3272 if (index(mapcard,'PHI').gt.0) then
3274 else if (index(mapcard,'THE').gt.0) then
3276 else if (index(mapcard,'ALP').gt.0) then
3278 else if (index(mapcard,'OME').gt.0) then
3281 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3282 stop 'Error - illegal variable spec in MAP card.'
3284 call readi (mapcard,'RES1',res1(imap),0)
3285 call readi (mapcard,'RES2',res2(imap),0)
3286 if (res1(imap).eq.0) then
3287 res1(imap)=res2(imap)
3288 else if (res2(imap).eq.0) then
3289 res2(imap)=res1(imap)
3291 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3292 write (iout,'(a)') &
3293 'Error - illegal definition of variable group in MAP.'
3294 stop 'Error - illegal definition of variable group in MAP.'
3296 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3297 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3298 call readi(mapcard,'NSTEP',nstep(imap),0)
3299 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3300 write (iout,'(a)') &
3301 'Illegal boundary and/or step size specification in MAP.'
3302 stop 'Illegal boundary and/or step size specification in MAP.'
3306 end subroutine map_read
3307 !-----------------------------------------------------------------------------
3310 use control_data, only: vdisulf
3312 ! implicit real*8 (a-h,o-z)
3313 ! include 'DIMENSIONS'
3314 ! include 'COMMON.IOUNITS'
3315 ! include 'COMMON.GEO'
3316 ! include 'COMMON.CSA'
3317 ! include 'COMMON.BANK'
3318 ! include 'COMMON.CONTROL'
3319 ! character(len=80) :: ucase
3320 character(len=620) :: mcmcard
3322 ! integer :: ntf,ik,iw_pdb
3323 ! real(kind=8) :: var,ene
3325 call card_concat(mcmcard,.true.)
3327 call readi(mcmcard,'NCONF',nconf,50)
3328 call readi(mcmcard,'NADD',nadd,0)
3329 call readi(mcmcard,'JSTART',jstart,1)
3330 call readi(mcmcard,'JEND',jend,1)
3331 call readi(mcmcard,'NSTMAX',nstmax,500000)
3332 call readi(mcmcard,'N0',n0,1)
3333 call readi(mcmcard,'N1',n1,6)
3334 call readi(mcmcard,'N2',n2,4)
3335 call readi(mcmcard,'N3',n3,0)
3336 call readi(mcmcard,'N4',n4,0)
3337 call readi(mcmcard,'N5',n5,0)
3338 call readi(mcmcard,'N6',n6,10)
3339 call readi(mcmcard,'N7',n7,0)
3340 call readi(mcmcard,'N8',n8,0)
3341 call readi(mcmcard,'N9',n9,0)
3342 call readi(mcmcard,'N14',n14,0)
3343 call readi(mcmcard,'N15',n15,0)
3344 call readi(mcmcard,'N16',n16,0)
3345 call readi(mcmcard,'N17',n17,0)
3346 call readi(mcmcard,'N18',n18,0)
3348 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3350 call readi(mcmcard,'NDIFF',ndiff,2)
3351 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3352 call readi(mcmcard,'IS1',is1,1)
3353 call readi(mcmcard,'IS2',is2,8)
3354 call readi(mcmcard,'NRAN0',nran0,4)
3355 call readi(mcmcard,'NRAN1',nran1,2)
3356 call readi(mcmcard,'IRR',irr,1)
3357 call readi(mcmcard,'NSEED',nseed,20)
3358 call readi(mcmcard,'NTOTAL',ntotal,10000)
3359 call reada(mcmcard,'CUT1',cut1,2.0d0)
3360 call reada(mcmcard,'CUT2',cut2,5.0d0)
3361 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3362 call readi(mcmcard,'ICMAX',icmax,3)
3363 call readi(mcmcard,'IRESTART',irestart,0)
3364 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3367 call reada(mcmcard,'DELE',dele,20.0d0)
3368 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3369 call readi(mcmcard,'IREF',iref,0)
3370 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3371 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3372 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3373 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3374 write (iout,*) "NCONF_IN",nconf_in
3376 end subroutine csaread
3377 !-----------------------------------------------------------------------------
3381 use control_data, only: MaxMoveType
3384 ! implicit real*8 (a-h,o-z)
3385 ! include 'DIMENSIONS'
3386 ! include 'COMMON.MCM'
3387 ! include 'COMMON.MCE'
3388 ! include 'COMMON.IOUNITS'
3389 ! character(len=80) :: ucase
3390 character(len=320) :: mcmcard
3393 ! real(kind=8) :: var,ene
3395 call card_concat(mcmcard,.true.)
3396 call readi(mcmcard,'MAXACC',maxacc,100)
3397 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3398 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3399 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3400 call readi(mcmcard,'MAXREPM',maxrepm,200)
3401 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3402 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3403 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3404 call reada(mcmcard,'E_UP',e_up,5.0D0)
3405 call reada(mcmcard,'DELTE',delte,0.1D0)
3406 call readi(mcmcard,'NSWEEP',nsweep,5)
3407 call readi(mcmcard,'NSTEPH',nsteph,0)
3408 call readi(mcmcard,'NSTEPC',nstepc,0)
3409 call reada(mcmcard,'TMIN',tmin,298.0D0)
3410 call reada(mcmcard,'TMAX',tmax,298.0D0)
3411 call readi(mcmcard,'NWINDOW',nwindow,0)
3412 call readi(mcmcard,'PRINT_MC',print_mc,0)
3413 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3414 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3415 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3416 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3417 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3418 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3419 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3420 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3421 if (nwindow.gt.0) then
3422 allocate(winstart(nwindow)) !!el (maxres)
3423 allocate(winend(nwindow)) !!el
3424 allocate(winlen(nwindow)) !!el
3425 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3427 winlen(i)=winend(i)-winstart(i)+1
3430 if (tmax.lt.tmin) tmax=tmin
3431 if (tmax.eq.tmin) then
3435 if (nstepc.gt.0 .and. nsteph.gt.0) then
3436 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3437 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3439 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3440 ! Probabilities of different move types
3441 sumpro_type(0)=0.0D0
3442 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3443 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3444 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3445 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3446 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3447 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3448 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3450 print *,'i',i,' sumprotype',sumpro_type(i)
3451 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3452 print *,'i',i,' sumprotype',sumpro_type(i)
3455 end subroutine mcmread
3456 !-----------------------------------------------------------------------------
3457 subroutine read_minim
3460 ! implicit real*8 (a-h,o-z)
3461 ! include 'DIMENSIONS'
3462 ! include 'COMMON.MINIM'
3463 ! include 'COMMON.IOUNITS'
3464 ! character(len=80) :: ucase
3465 character(len=320) :: minimcard
3467 ! integer :: ntf,ik,iw_pdb
3468 ! real(kind=8) :: var,ene
3470 call card_concat(minimcard,.true.)
3471 call readi(minimcard,'MAXMIN',maxmin,2000)
3472 call readi(minimcard,'MAXFUN',maxfun,5000)
3473 call readi(minimcard,'MINMIN',minmin,maxmin)
3474 call readi(minimcard,'MINFUN',minfun,maxmin)
3475 call reada(minimcard,'TOLF',tolf,1.0D-2)
3476 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3477 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3478 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3479 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3480 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3481 'Options in energy minimization:'
3482 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3483 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3484 'MinMin:',MinMin,' MinFun:',MinFun,&
3485 ' TolF:',TolF,' RTolF:',RTolF
3487 end subroutine read_minim
3488 !-----------------------------------------------------------------------------
3489 subroutine openunits
3491 use MD_data, only: usampl
3494 use control_data, only:out1file
3495 use control, only: getenv_loc
3496 ! implicit real*8 (a-h,o-z)
3497 ! include 'DIMENSIONS'
3500 character(len=16) :: form,nodename
3501 integer :: nodelen,ierror,npos
3503 ! include 'COMMON.SETUP'
3504 ! include 'COMMON.IOUNITS'
3505 ! include 'COMMON.MD'
3506 ! include 'COMMON.CONTROL'
3507 integer :: lenpre,lenpot,lentmp !,ilen
3509 character(len=3) :: out1file_text !,ucase
3510 character(len=3) :: ll
3513 ! integer :: ntf,ik,iw_pdb
3514 ! real(kind=8) :: var,ene
3516 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3517 call getenv_loc("PREFIX",prefix)
3519 call getenv_loc("POT",pot)
3520 call getenv_loc("DIRTMP",tmpdir)
3521 call getenv_loc("CURDIR",curdir)
3522 call getenv_loc("OUT1FILE",out1file_text)
3523 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3524 out1file_text=ucase(out1file_text)
3525 if (out1file_text(1:1).eq."Y") then
3528 out1file=fg_rank.gt.0
3533 if (lentmp.gt.0) then
3534 write (*,'(80(1h!))')
3535 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3536 write (*,'(80(1h!))')
3537 write (*,*)"All output files will be on node /tmp directory."
3539 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3540 if (me.eq.king) then
3541 write (*,*) "The master node is ",nodename
3542 else if (fg_rank.eq.0) then
3543 write (*,*) "I am the CG slave node ",nodename
3545 write (*,*) "I am the FG slave node ",nodename
3548 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3549 lenpre = lentmp+lenpre+1
3551 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3552 ! Get the names and open the input files
3553 #if defined(WINIFL) || defined(WINPGI)
3554 open(1,file=pref_orig(:ilen(pref_orig))// &
3555 '.inp',status='old',readonly,shared)
3556 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3557 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3558 ! Get parameter filenames and open the parameter files.
3559 call getenv_loc('BONDPAR',bondname)
3560 open (ibond,file=bondname,status='old',readonly,shared)
3561 call getenv_loc('THETPAR',thetname)
3562 open (ithep,file=thetname,status='old',readonly,shared)
3563 call getenv_loc('ROTPAR',rotname)
3564 open (irotam,file=rotname,status='old',readonly,shared)
3565 call getenv_loc('TORPAR',torname)
3566 open (itorp,file=torname,status='old',readonly,shared)
3567 call getenv_loc('TORDPAR',tordname)
3568 open (itordp,file=tordname,status='old',readonly,shared)
3569 call getenv_loc('FOURIER',fouriername)
3570 open (ifourier,file=fouriername,status='old',readonly,shared)
3571 call getenv_loc('ELEPAR',elename)
3572 open (ielep,file=elename,status='old',readonly,shared)
3573 call getenv_loc('SIDEPAR',sidename)
3574 open (isidep,file=sidename,status='old',readonly,shared)
3575 #elif (defined CRAY) || (defined AIX)
3576 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3578 ! print *,"Processor",myrank," opened file 1"
3579 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3580 ! print *,"Processor",myrank," opened file 9"
3581 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3582 ! Get parameter filenames and open the parameter files.
3583 call getenv_loc('BONDPAR',bondname)
3584 open (ibond,file=bondname,status='old',action='read')
3585 ! print *,"Processor",myrank," opened file IBOND"
3586 call getenv_loc('THETPAR',thetname)
3587 open (ithep,file=thetname,status='old',action='read')
3588 ! print *,"Processor",myrank," opened file ITHEP"
3589 call getenv_loc('ROTPAR',rotname)
3590 open (irotam,file=rotname,status='old',action='read')
3591 ! print *,"Processor",myrank," opened file IROTAM"
3592 call getenv_loc('TORPAR',torname)
3593 open (itorp,file=torname,status='old',action='read')
3594 ! print *,"Processor",myrank," opened file ITORP"
3595 call getenv_loc('TORDPAR',tordname)
3596 open (itordp,file=tordname,status='old',action='read')
3597 ! print *,"Processor",myrank," opened file ITORDP"
3598 call getenv_loc('SCCORPAR',sccorname)
3599 open (isccor,file=sccorname,status='old',action='read')
3600 ! print *,"Processor",myrank," opened file ISCCOR"
3601 call getenv_loc('FOURIER',fouriername)
3602 open (ifourier,file=fouriername,status='old',action='read')
3603 ! print *,"Processor",myrank," opened file IFOURIER"
3604 call getenv_loc('ELEPAR',elename)
3605 open (ielep,file=elename,status='old',action='read')
3606 ! print *,"Processor",myrank," opened file IELEP"
3607 call getenv_loc('SIDEPAR',sidename)
3608 open (isidep,file=sidename,status='old',action='read')
3609 ! print *,"Processor",myrank," opened file ISIDEP"
3610 ! print *,"Processor",myrank," opened parameter files"
3612 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3613 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3614 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3615 ! Get parameter filenames and open the parameter files.
3616 call getenv_loc('BONDPAR',bondname)
3617 open (ibond,file=bondname,status='old')
3618 call getenv_loc('THETPAR',thetname)
3619 open (ithep,file=thetname,status='old')
3620 call getenv_loc('ROTPAR',rotname)
3621 open (irotam,file=rotname,status='old')
3622 call getenv_loc('TORPAR',torname)
3623 open (itorp,file=torname,status='old')
3624 call getenv_loc('TORDPAR',tordname)
3625 open (itordp,file=tordname,status='old')
3626 call getenv_loc('SCCORPAR',sccorname)
3627 open (isccor,file=sccorname,status='old')
3628 call getenv_loc('FOURIER',fouriername)
3629 open (ifourier,file=fouriername,status='old')
3630 call getenv_loc('ELEPAR',elename)
3631 open (ielep,file=elename,status='old')
3632 call getenv_loc('SIDEPAR',sidename)
3633 open (isidep,file=sidename,status='old')
3635 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3637 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3638 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3639 ! Get parameter filenames and open the parameter files.
3640 call getenv_loc('BONDPAR',bondname)
3641 open (ibond,file=bondname,status='old',action='read')
3642 call getenv_loc('THETPAR',thetname)
3643 open (ithep,file=thetname,status='old',action='read')
3644 call getenv_loc('ROTPAR',rotname)
3645 open (irotam,file=rotname,status='old',action='read')
3646 call getenv_loc('TORPAR',torname)
3647 open (itorp,file=torname,status='old',action='read')
3648 call getenv_loc('TORDPAR',tordname)
3649 open (itordp,file=tordname,status='old',action='read')
3650 call getenv_loc('SCCORPAR',sccorname)
3651 open (isccor,file=sccorname,status='old',action='read')
3653 call getenv_loc('THETPARPDB',thetname_pdb)
3654 print *,"thetname_pdb ",thetname_pdb
3655 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3656 print *,ithep_pdb," opened"
3658 call getenv_loc('FOURIER',fouriername)
3659 open (ifourier,file=fouriername,status='old',readonly)
3660 call getenv_loc('ELEPAR',elename)
3661 open (ielep,file=elename,status='old',readonly)
3662 call getenv_loc('SIDEPAR',sidename)
3663 open (isidep,file=sidename,status='old',readonly)
3665 call getenv_loc('ROTPARPDB',rotname_pdb)
3666 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3671 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3672 ! Use -DOLDSCP to use hard-coded constants instead.
3674 call getenv_loc('SCPPAR',scpname)
3675 #if defined(WINIFL) || defined(WINPGI)
3676 open (iscpp,file=scpname,status='old',readonly,shared)
3677 #elif (defined CRAY) || (defined AIX)
3678 open (iscpp,file=scpname,status='old',action='read')
3680 open (iscpp,file=scpname,status='old')
3682 open (iscpp,file=scpname,status='old',action='read')
3685 call getenv_loc('PATTERN',patname)
3686 #if defined(WINIFL) || defined(WINPGI)
3687 open (icbase,file=patname,status='old',readonly,shared)
3688 #elif (defined CRAY) || (defined AIX)
3689 open (icbase,file=patname,status='old',action='read')
3691 open (icbase,file=patname,status='old')
3693 open (icbase,file=patname,status='old',action='read')
3696 ! Open output file only for CG processes
3697 ! print *,"Processor",myrank," fg_rank",fg_rank
3698 if (fg_rank.eq.0) then
3700 if (nodes.eq.1) then
3703 npos = dlog10(dfloat(nodes-1))+1
3705 if (npos.lt.3) npos=3
3706 write (liczba,'(i1)') npos
3707 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3709 write (liczba,form) me
3710 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3711 liczba(:ilen(liczba))
3712 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3714 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3716 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3717 liczba(:ilen(liczba))//'.mol2'
3718 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3719 liczba(:ilen(liczba))//'.stat'
3721 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3722 //liczba(:ilen(liczba))//'.stat')
3723 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3726 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3727 liczba(:ilen(liczba))//'.const'
3732 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3733 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3734 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3735 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3736 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3738 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3740 rest2name=prefix(:ilen(prefix))//'.rst'
3742 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3745 #if defined(AIX) || defined(PGI)
3746 if (me.eq.king .or. .not. out1file) &
3747 open(iout,file=outname,status='unknown')
3749 if (fg_rank.gt.0) then
3750 write (liczba,'(i3.3)') myrank/nfgtasks
3751 write (ll,'(bz,i3.3)') fg_rank
3752 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3757 open(igeom,file=intname,status='unknown',position='append')
3758 open(ipdb,file=pdbname,status='unknown')
3759 open(imol2,file=mol2name,status='unknown')
3760 open(istat,file=statname,status='unknown',position='append')
3762 !1out open(iout,file=outname,status='unknown')
3765 if (me.eq.king .or. .not.out1file) &
3766 open(iout,file=outname,status='unknown')
3768 if (fg_rank.gt.0) then
3769 write (liczba,'(i3.3)') myrank/nfgtasks
3770 write (ll,'(bz,i3.3)') fg_rank
3771 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3776 open(igeom,file=intname,status='unknown',access='append')
3777 open(ipdb,file=pdbname,status='unknown')
3778 open(imol2,file=mol2name,status='unknown')
3779 open(istat,file=statname,status='unknown',access='append')
3781 !1out open(iout,file=outname,status='unknown')
3784 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3785 csa_seed=prefix(:lenpre)//'.CSA.seed'
3786 csa_history=prefix(:lenpre)//'.CSA.history'
3787 csa_bank=prefix(:lenpre)//'.CSA.bank'
3788 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3789 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3790 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3791 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3792 csa_int=prefix(:lenpre)//'.int'
3793 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3794 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3795 csa_in=prefix(:lenpre)//'.CSA.in'
3796 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3799 write (iout,'(80(1h-))')
3800 write (iout,'(30x,a)') "FILE ASSIGNMENT"
3801 write (iout,'(80(1h-))')
3802 write (iout,*) "Input file : ",&
3803 pref_orig(:ilen(pref_orig))//'.inp'
3804 write (iout,*) "Output file : ",&
3805 outname(:ilen(outname))
3807 write (iout,*) "Sidechain potential file : ",&
3808 sidename(:ilen(sidename))
3810 write (iout,*) "SCp potential file : ",&
3811 scpname(:ilen(scpname))
3813 write (iout,*) "Electrostatic potential file : ",&
3814 elename(:ilen(elename))
3815 write (iout,*) "Cumulant coefficient file : ",&
3816 fouriername(:ilen(fouriername))
3817 write (iout,*) "Torsional parameter file : ",&
3818 torname(:ilen(torname))
3819 write (iout,*) "Double torsional parameter file : ",&
3820 tordname(:ilen(tordname))
3821 write (iout,*) "SCCOR parameter file : ",&
3822 sccorname(:ilen(sccorname))
3823 write (iout,*) "Bond & inertia constant file : ",&
3824 bondname(:ilen(bondname))
3825 write (iout,*) "Bending parameter file : ",&
3826 thetname(:ilen(thetname))
3827 write (iout,*) "Rotamer parameter file : ",&
3828 rotname(:ilen(rotname))
3831 write (iout,*) "Thetpdb parameter file : ",&
3832 thetname_pdb(:ilen(thetname_pdb))
3835 write (iout,*) "Threading database : ",&
3836 patname(:ilen(patname))
3838 write (iout,*)" DIRTMP : ",&
3840 write (iout,'(80(1h-))')
3843 end subroutine openunits
3844 !-----------------------------------------------------------------------------
3847 use geometry_data, only: nres,dc
3849 ! implicit real*8 (a-h,o-z)
3850 ! include 'DIMENSIONS'
3851 ! include 'COMMON.CHAIN'
3852 ! include 'COMMON.IOUNITS'
3853 ! include 'COMMON.MD'
3856 ! real(kind=8) :: var,ene
3858 open(irest2,file=rest2name,status='unknown')
3859 read(irest2,*) totT,EK,potE,totE,t_bath
3861 ! AL 4/17/17: Now reading d_t(0,:) too
3863 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
3866 ! AL 4/17/17: Now reading d_c(0,:) too
3868 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
3871 read (irest2,*) iset
3875 end subroutine readrst
3876 !-----------------------------------------------------------------------------
3877 subroutine read_fragments
3881 use control_data, only:out1file
3884 ! implicit real*8 (a-h,o-z)
3885 ! include 'DIMENSIONS'
3889 ! include 'COMMON.SETUP'
3890 ! include 'COMMON.CHAIN'
3891 ! include 'COMMON.IOUNITS'
3892 ! include 'COMMON.MD'
3893 ! include 'COMMON.CONTROL'
3896 ! real(kind=8) :: var,ene
3898 read(inp,*) nset,nfrag,npair,nfrag_back
3900 !el from module energy
3901 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
3902 if(.not.allocated(wfrag_back)) then
3903 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3904 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3906 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
3907 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
3909 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
3910 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
3913 if(me.eq.king.or..not.out1file) &
3914 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
3915 " nfrag_back",nfrag_back
3917 read(inp,*) mset(iset)
3919 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
3921 if(me.eq.king.or..not.out1file) &
3922 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
3923 ifrag(2,i,iset), qinfrag(i,iset)
3926 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
3928 if(me.eq.king.or..not.out1file) &
3929 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
3930 ipair(2,i,iset), qinpair(i,iset)
3933 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3934 wfrag_back(3,i,iset),&
3935 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
3936 if(me.eq.king.or..not.out1file) &
3937 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3938 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
3942 end subroutine read_fragments
3943 !-----------------------------------------------------------------------------
3945 !-----------------------------------------------------------------------------
3949 ! implicit real*8 (a-h,o-z)
3950 ! include 'DIMENSIONS'
3951 ! include 'COMMON.CSA'
3952 ! include 'COMMON.BANK'
3953 ! include 'COMMON.IOUNITS'
3955 ! integer :: ntf,ik,iw_pdb
3956 ! real(kind=8) :: var,ene
3958 open(icsa_in,file=csa_in,status="old",err=100)
3959 read(icsa_in,*) nconf
3960 read(icsa_in,*) jstart,jend
3961 read(icsa_in,*) nstmax
3962 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
3963 read(icsa_in,*) nran0,nran1,irr
3964 read(icsa_in,*) nseed
3965 read(icsa_in,*) ntotal,cut1,cut2
3966 read(icsa_in,*) estop
3967 read(icsa_in,*) icmax,irestart
3968 read(icsa_in,*) ntbankm,dele,difcut
3969 read(icsa_in,*) iref,rmscut,pnccut
3970 read(icsa_in,*) ndiff
3977 end subroutine csa_read
3978 !-----------------------------------------------------------------------------
3979 subroutine initial_write
3982 ! implicit real*8 (a-h,o-z)
3983 ! include 'DIMENSIONS'
3984 ! include 'COMMON.CSA'
3985 ! include 'COMMON.BANK'
3986 ! include 'COMMON.IOUNITS'
3988 ! integer :: ntf,ik,iw_pdb
3989 ! real(kind=8) :: var,ene
3991 open(icsa_seed,file=csa_seed,status="unknown")
3992 write(icsa_seed,*) "seed"
3994 #if defined(AIX) || defined(PGI)
3995 open(icsa_history,file=csa_history,status="unknown",&
3998 open(icsa_history,file=csa_history,status="unknown",&
4001 write(icsa_history,*) nconf
4002 write(icsa_history,*) jstart,jend
4003 write(icsa_history,*) nstmax
4004 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4005 write(icsa_history,*) nran0,nran1,irr
4006 write(icsa_history,*) nseed
4007 write(icsa_history,*) ntotal,cut1,cut2
4008 write(icsa_history,*) estop
4009 write(icsa_history,*) icmax,irestart
4010 write(icsa_history,*) ntbankm,dele,difcut
4011 write(icsa_history,*) iref,rmscut,pnccut
4012 write(icsa_history,*) ndiff
4014 write(icsa_history,*)
4017 open(icsa_bank1,file=csa_bank1,status="unknown")
4018 write(icsa_bank1,*) 0
4022 end subroutine initial_write
4023 !-----------------------------------------------------------------------------
4024 subroutine restart_write
4027 ! implicit real*8 (a-h,o-z)
4028 ! include 'DIMENSIONS'
4029 ! include 'COMMON.IOUNITS'
4030 ! include 'COMMON.CSA'
4031 ! include 'COMMON.BANK'
4033 ! integer :: ntf,ik,iw_pdb
4034 ! real(kind=8) :: var,ene
4036 #if defined(AIX) || defined(PGI)
4037 open(icsa_history,file=csa_history,position="append")
4039 open(icsa_history,file=csa_history,access="append")
4041 write(icsa_history,*)
4042 write(icsa_history,*) "This is restart"
4043 write(icsa_history,*)
4044 write(icsa_history,*) nconf
4045 write(icsa_history,*) jstart,jend
4046 write(icsa_history,*) nstmax
4047 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4048 write(icsa_history,*) nran0,nran1,irr
4049 write(icsa_history,*) nseed
4050 write(icsa_history,*) ntotal,cut1,cut2
4051 write(icsa_history,*) estop
4052 write(icsa_history,*) icmax,irestart
4053 write(icsa_history,*) ntbankm,dele,difcut
4054 write(icsa_history,*) iref,rmscut,pnccut
4055 write(icsa_history,*) ndiff
4056 write(icsa_history,*)
4057 write(icsa_history,*) "irestart is: ", irestart
4059 write(icsa_history,*)
4063 end subroutine restart_write
4064 !-----------------------------------------------------------------------------
4066 !-----------------------------------------------------------------------------
4067 subroutine write_pdb(npdb,titelloc,ee)
4069 ! implicit real*8 (a-h,o-z)
4070 ! include 'DIMENSIONS'
4071 ! include 'COMMON.IOUNITS'
4072 character(len=50) :: titelloc1
4073 character*(*) :: titelloc
4074 character(len=3) :: zahl
4075 character(len=5) :: liczba5
4077 integer :: npdb !,ilen
4081 ! real(kind=8) :: var,ene
4085 if (npdb.lt.1000) then
4086 call numstr(npdb,zahl)
4087 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4089 if (npdb.lt.10000) then
4090 write(liczba5,'(i1,i4)') 0,npdb
4092 write(liczba5,'(i5)') npdb
4094 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4096 call pdbout(ee,titelloc1,ipdb)
4099 end subroutine write_pdb
4100 !-----------------------------------------------------------------------------
4102 !-----------------------------------------------------------------------------
4103 subroutine write_thread_summary
4104 ! Thread the sequence through a database of known structures
4105 use control_data, only: refstr
4107 use energy_data, only: n_ene_comp
4109 ! implicit real*8 (a-h,o-z)
4110 ! include 'DIMENSIONS'
4112 use MPI_data !include 'COMMON.INFO'
4115 ! include 'COMMON.CONTROL'
4116 ! include 'COMMON.CHAIN'
4117 ! include 'COMMON.DBASE'
4118 ! include 'COMMON.INTERACT'
4119 ! include 'COMMON.VAR'
4120 ! include 'COMMON.THREAD'
4121 ! include 'COMMON.FFIELD'
4122 ! include 'COMMON.SBRIDGE'
4123 ! include 'COMMON.HEADER'
4124 ! include 'COMMON.NAMES'
4125 ! include 'COMMON.IOUNITS'
4126 ! include 'COMMON.TIME1'
4128 integer,dimension(maxthread) :: ip
4129 real(kind=8),dimension(0:n_ene) :: energia
4131 integer :: i,j,ii,jj,ipj,ik,kk,ist
4132 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4134 write (iout,'(30x,a/)') &
4135 ' *********** Summary threading statistics ************'
4136 write (iout,'(a)') 'Initial energies:'
4137 write (iout,'(a4,2x,a12,14a14,3a8)') &
4138 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4139 'RMSnat','NatCONT','NNCONT','RMS'
4140 ! Energy sort patterns
4145 enet=ener(n_ene-1,ip(i))
4148 if (ener(n_ene-1,ip(j)).lt.enet) then
4150 enet=ener(n_ene-1,ip(j))
4162 ist=nres_base(2,ii)+ipatt(2,i)
4164 energia(i)=ener0(kk,i)
4166 etot=ener0(n_ene_comp+1,i)
4167 rmsnat=ener0(n_ene_comp+2,i)
4168 rms=ener0(n_ene_comp+3,i)
4169 frac=ener0(n_ene_comp+4,i)
4170 frac_nn=ener0(n_ene_comp+5,i)
4173 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4174 i,str_nam(ii),ist+1,&
4175 (energia(print_order(kk)),kk=1,nprint_ene),&
4176 etot,rmsnat,frac,frac_nn,rms
4178 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4179 i,str_nam(ii),ist+1,&
4180 (energia(print_order(kk)),kk=1,nprint_ene),etot
4183 write (iout,'(//a)') 'Final energies:'
4184 write (iout,'(a4,2x,a12,17a14,3a8)') &
4185 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4186 'RMSnat','NatCONT','NNCONT','RMS'
4190 ist=nres_base(2,ii)+ipatt(2,i)
4192 energia(kk)=ener(kk,ik)
4194 etot=ener(n_ene_comp+1,i)
4195 rmsnat=ener(n_ene_comp+2,i)
4196 rms=ener(n_ene_comp+3,i)
4197 frac=ener(n_ene_comp+4,i)
4198 frac_nn=ener(n_ene_comp+5,i)
4199 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4200 i,str_nam(ii),ist+1,&
4201 (energia(print_order(kk)),kk=1,nprint_ene),&
4202 etot,rmsnat,frac,frac_nn,rms
4204 write (iout,'(/a/)') 'IEXAM array:'
4205 write (iout,'(i5)') nexcl
4207 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4209 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4210 'Max. time for threading step ',max_time_for_thread,&
4211 'Average time for threading step: ',ave_time_for_thread
4213 end subroutine write_thread_summary
4214 !-----------------------------------------------------------------------------
4215 subroutine write_stat_thread(ithread,ipattern,ist)
4217 use energy_data, only: n_ene_comp
4219 ! implicit real*8 (a-h,o-z)
4220 ! include "DIMENSIONS"
4221 ! include "COMMON.CONTROL"
4222 ! include "COMMON.IOUNITS"
4223 ! include "COMMON.THREAD"
4224 ! include "COMMON.FFIELD"
4225 ! include "COMMON.DBASE"
4226 ! include "COMMON.NAMES"
4227 real(kind=8),dimension(0:n_ene) :: energia
4229 integer :: ithread,ipattern,ist,i
4230 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4232 #if defined(AIX) || defined(PGI)
4233 open(istat,file=statname,position='append')
4235 open(istat,file=statname,access='append')
4238 energia(i)=ener(i,ithread)
4240 etot=ener(n_ene_comp+1,ithread)
4241 rmsnat=ener(n_ene_comp+2,ithread)
4242 rms=ener(n_ene_comp+3,ithread)
4243 frac=ener(n_ene_comp+4,ithread)
4244 frac_nn=ener(n_ene_comp+5,ithread)
4245 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4246 ithread,str_nam(ipattern),ist+1,&
4247 (energia(print_order(i)),i=1,nprint_ene),&
4248 etot,rmsnat,frac,frac_nn,rms
4251 end subroutine write_stat_thread
4252 !-----------------------------------------------------------------------------
4254 !-----------------------------------------------------------------------------
4255 end module io_config