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:maxtor,maxterm
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)
840 athet(j,i,ichir1,ichir2)=0.0D0
841 bthet(j,i,ichir1,ichir2)=0.0D0
858 ! Read the parameters of the probability distribution/energy expression
859 ! of the virtual-bond valence angles theta
862 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
863 (bthet(j,i,1,1),j=1,2)
864 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
865 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
866 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
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)
881 athet(1,i,-1,-1)=athet(1,-i,1,1)
882 athet(2,i,-1,-1)=-athet(2,-i,1,1)
883 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
884 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
885 athet(1,i,-1,1)=athet(1,-i,1,1)
886 athet(2,i,-1,1)=-athet(2,-i,1,1)
887 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
888 bthet(2,i,-1,1)=bthet(2,-i,1,1)
889 athet(1,i,1,-1)=-athet(1,-i,1,1)
890 athet(2,i,1,-1)=athet(2,-i,1,1)
891 bthet(1,i,1,-1)=bthet(1,-i,1,1)
892 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
897 polthet(j,i)=polthet(j,-i)
900 gthet(j,i)=gthet(j,-i)
908 'Parameters of the virtual-bond valence angles:'
909 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
910 ' ATHETA0 ',' A1 ',' A2 ',&
913 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
914 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
916 write (iout,'(/a/9x,5a/79(1h-))') &
917 'Parameters of the expression for sigma(theta_c):',&
918 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
919 ' ALPH3 ',' SIGMA0C '
921 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
922 (polthet(j,i),j=0,3),sigc0(i)
924 write (iout,'(/a/9x,5a/79(1h-))') &
925 'Parameters of the second gaussian:',&
926 ' THETA0 ',' SIGMA0 ',' G1 ',&
929 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
930 sig0(i),(gthet(j,i),j=1,3)
934 'Parameters of the virtual-bond valence angles:'
935 write (iout,'(/a/9x,5a/79(1h-))') &
936 'Coefficients of expansion',&
937 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
938 ' b1*10^1 ',' b2*10^1 '
940 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
941 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
942 (10*bthet(j,i,1,1),j=1,2)
944 write (iout,'(/a/9x,5a/79(1h-))') &
945 'Parameters of the expression for sigma(theta_c):',&
946 ' alpha0 ',' alph1 ',' alph2 ',&
947 ' alhp3 ',' sigma0c '
949 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
950 (polthet(j,i),j=0,3),sigc0(i)
952 write (iout,'(/a/9x,5a/79(1h-))') &
953 'Parameters of the second gaussian:',&
954 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
957 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
958 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
964 ! Read the parameters of Utheta determined from ab initio surfaces
965 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
967 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
968 ntheterm3,nsingle,ndouble
969 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
971 !----------------------------------------------------
972 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
973 allocate(aa0thet(-maxthetyp1:maxthetyp1,&
974 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
975 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
976 allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
977 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
978 !(maxtheterm,-maxthetyp1:maxthetyp1,&
979 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
980 allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
981 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
982 allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
983 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
984 allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
985 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
986 allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
987 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
988 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
989 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
990 allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
991 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
992 allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
993 -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
994 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
995 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
997 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
999 ithetyp(i)=-ithetyp(-i)
1002 do i=-maxthetyp,maxthetyp
1003 do j=-maxthetyp,maxthetyp
1004 do k=-maxthetyp,maxthetyp
1005 aa0thet(i,j,k,iblock)=0.0d0
1007 aathet(l,i,j,k,iblock)=0.0d0
1011 bbthet(m,l,i,j,k,iblock)=0.0d0
1012 ccthet(m,l,i,j,k,iblock)=0.0d0
1013 ddthet(m,l,i,j,k,iblock)=0.0d0
1014 eethet(m,l,i,j,k,iblock)=0.0d0
1020 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
1021 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
1029 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1031 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1032 ! VAR:1=non-glicyne non-proline 2=proline
1033 ! VAR:negative values for D-aminoacid
1035 do j=-nthetyp,nthetyp
1036 do k=-nthetyp,nthetyp
1037 read (ithep,'(6a)',end=111,err=111) res1
1038 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1039 ! VAR: aa0thet is variable describing the average value of Foureir
1040 ! VAR: expansion series
1041 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1042 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1043 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1044 read (ithep,*,end=111,err=111) &
1045 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1046 read (ithep,*,end=111,err=111) &
1047 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1048 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1049 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1050 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1052 read (ithep,*,end=111,err=111) &
1053 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1054 ffthet(lll,llll,ll,i,j,k,iblock),&
1055 ggthet(llll,lll,ll,i,j,k,iblock),&
1056 ggthet(lll,llll,ll,i,j,k,iblock),&
1057 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1062 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1063 ! coefficients of theta-and-gamma-dependent terms are zero.
1064 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1065 ! RECOMENTDED AFTER VERSION 3.3)
1069 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1070 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1072 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1073 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1076 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1078 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1081 ! AND COMMENT THE LOOPS BELOW
1085 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1086 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1088 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1089 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1092 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1094 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1099 ! Substitution for D aminoacids from symmetry.
1102 do j=-nthetyp,nthetyp
1103 do k=-nthetyp,nthetyp
1104 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1106 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1110 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1111 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1112 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1113 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1119 ffthet(llll,lll,ll,i,j,k,iblock)= &
1120 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1121 ffthet(lll,llll,ll,i,j,k,iblock)= &
1122 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1123 ggthet(llll,lll,ll,i,j,k,iblock)= &
1124 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1125 ggthet(lll,llll,ll,i,j,k,iblock)= &
1126 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1135 ! Control printout of the coefficients of virtual-bond-angle potentials
1138 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1143 write (iout,'(//4a)') &
1144 'Type ',onelett(i),onelett(j),onelett(k)
1145 write (iout,'(//a,10x,a)') " l","a[l]"
1146 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1147 write (iout,'(i2,1pe15.5)') &
1148 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1150 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1151 "b",l,"c",l,"d",l,"e",l
1153 write (iout,'(i2,4(1pe15.5))') m,&
1154 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1155 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1159 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1160 "f+",l,"f-",l,"g+",l,"g-",l
1163 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1164 ffthet(n,m,l,i,j,k,iblock),&
1165 ffthet(m,n,l,i,j,k,iblock),&
1166 ggthet(n,m,l,i,j,k,iblock),&
1167 ggthet(m,n,l,i,j,k,iblock)
1177 write (2,*) "Start reading THETA_PDB",ithep_pdb
1179 ! write (2,*) 'i=',i
1180 read (ithep_pdb,*,err=111,end=111) &
1181 a0thet(i),(athet(j,i,1,1),j=1,2),&
1182 (bthet(j,i,1,1),j=1,2)
1183 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1184 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1185 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1186 sigc0(i)=sigc0(i)**2
1189 athet(1,i,1,-1)=athet(1,i,1,1)
1190 athet(2,i,1,-1)=athet(2,i,1,1)
1191 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1192 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1193 athet(1,i,-1,1)=-athet(1,i,1,1)
1194 athet(2,i,-1,1)=-athet(2,i,1,1)
1195 bthet(1,i,-1,1)=bthet(1,i,1,1)
1196 bthet(2,i,-1,1)=bthet(2,i,1,1)
1199 a0thet(i)=a0thet(-i)
1200 athet(1,i,-1,-1)=athet(1,-i,1,1)
1201 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1202 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1203 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1204 athet(1,i,-1,1)=athet(1,-i,1,1)
1205 athet(2,i,-1,1)=-athet(2,-i,1,1)
1206 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1207 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1208 athet(1,i,1,-1)=-athet(1,-i,1,1)
1209 athet(2,i,1,-1)=athet(2,-i,1,1)
1210 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1211 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1212 theta0(i)=theta0(-i)
1216 polthet(j,i)=polthet(j,-i)
1219 gthet(j,i)=gthet(j,-i)
1222 write (2,*) "End reading THETA_PDB"
1227 !-------------------------------------------
1228 allocate(nlob(ntyp1)) !(ntyp1)
1229 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1230 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1231 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1249 gaussc(l,k,j,i)=0.0D0
1257 ! Read the parameters of the probability distribution/energy expression
1258 ! of the side chains.
1261 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1265 dsc_inv(i)=1.0D0/dsc(i)
1276 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1277 ((blower(k,l,1),l=1,k),k=1,3)
1278 censc(1,1,-i)=censc(1,1,i)
1279 censc(2,1,-i)=censc(2,1,i)
1280 censc(3,1,-i)=-censc(3,1,i)
1282 read (irotam,*,end=112,err=112) bsc(j,i)
1283 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1284 ((blower(k,l,j),l=1,k),k=1,3)
1285 censc(1,j,-i)=censc(1,j,i)
1286 censc(2,j,-i)=censc(2,j,i)
1287 censc(3,j,-i)=-censc(3,j,i)
1288 ! BSC is amplitude of Gaussian
1295 akl=akl+blower(k,m,j)*blower(l,m,j)
1299 if (((k.eq.3).and.(l.ne.3)) &
1300 .or.((l.eq.3).and.(k.ne.3))) then
1301 gaussc(k,l,j,-i)=-akl
1302 gaussc(l,k,j,-i)=-akl
1304 gaussc(k,l,j,-i)=akl
1305 gaussc(l,k,j,-i)=akl
1314 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1317 if (nlobi.gt.0) then
1319 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1320 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1321 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1322 'log h',(bsc(j,i),j=1,nlobi)
1323 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1324 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1326 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1327 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1330 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1331 write (iout,'(a,f10.4,4(16x,f10.4))') &
1332 'Center ',(bsc(j,i),j=1,nlobi)
1333 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1342 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1343 ! added by Urszula Kozlowska 07/11/2007
1345 !el Maximum number of SC local term fitting function coefficiants
1346 !el integer,parameter :: maxsccoef=65
1348 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1351 read (irotam,*,end=112,err=112)
1353 read (irotam,*,end=112,err=112)
1356 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1361 ! Read the parameters of the probability distribution/energy expression
1362 ! of the side chains.
1364 write (2,*) "Start reading ROTAM_PDB"
1366 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1370 dsc_inv(i)=1.0D0/dsc(i)
1381 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1382 ((blower(k,l,1),l=1,k),k=1,3)
1384 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1385 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1386 ((blower(k,l,j),l=1,k),k=1,3)
1393 akl=akl+blower(k,m,j)*blower(l,m,j)
1403 write (2,*) "End reading ROTAM_PDB"
1409 ! Read torsional parameters in old format
1411 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1413 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1414 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1415 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1417 !el from energy module--------
1418 allocate(v1(nterm_old,ntortyp,ntortyp))
1419 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1420 !el---------------------------
1425 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1431 write (iout,'(/a/)') 'Torsional constants:'
1434 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1435 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1441 ! Read torsional parameters
1443 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1444 read (itorp,*,end=113,err=113) ntortyp
1445 !el from energy module---------
1446 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1447 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1449 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1450 allocate(vlor2(maxlor,ntortyp,ntortyp))
1451 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1452 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1454 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1455 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1456 !el---------------------------
1458 do i=-ntortyp,ntortyp
1459 do j=-ntortyp,ntortyp
1465 !el---------------------------
1467 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1469 itortyp(i)=-itortyp(-i)
1471 ! itortyp(ntyp1)=ntortyp
1472 ! itortyp(-ntyp1)=-ntortyp
1474 write (iout,*) 'ntortyp',ntortyp
1476 do j=-ntortyp+1,ntortyp-1
1477 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1479 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1480 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1483 do k=1,nterm(i,j,iblock)
1484 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1486 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1487 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1488 v0ij=v0ij+si*v1(k,i,j,iblock)
1490 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1491 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1492 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1494 do k=1,nlor(i,j,iblock)
1495 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1496 vlor2(k,i,j),vlor3(k,i,j)
1497 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1500 v0(-i,-j,iblock)=v0ij
1506 write (iout,'(/a/)') 'Torsional constants:'
1508 do i=-ntortyp,ntortyp
1509 do j=-ntortyp,ntortyp
1510 write (iout,*) 'ityp',i,' jtyp',j
1511 write (iout,*) 'Fourier constants'
1512 do k=1,nterm(i,j,iblock)
1513 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1516 write (iout,*) 'Lorenz constants'
1517 do k=1,nlor(i,j,iblock)
1518 write (iout,'(3(1pe15.5))') &
1519 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1525 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1527 ! 6/23/01 Read parameters for double torsionals
1529 !el from energy module------------
1530 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1531 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1532 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1533 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1534 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1535 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1536 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1537 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1538 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1539 !---------------------------------
1543 do j=-ntortyp+1,ntortyp-1
1544 do k=-ntortyp+1,ntortyp-1
1545 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1546 ! write (iout,*) "OK onelett",
1549 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1550 .or. t3.ne.toronelet(k)) then
1551 write (iout,*) "Error in double torsional parameter file",&
1554 call MPI_Finalize(Ierror)
1556 stop "Error in double torsional parameter file"
1558 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1559 ntermd_2(i,j,k,iblock)
1560 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1561 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1562 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1563 ntermd_1(i,j,k,iblock))
1564 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1565 ntermd_1(i,j,k,iblock))
1566 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1567 ntermd_1(i,j,k,iblock))
1568 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1569 ntermd_1(i,j,k,iblock))
1570 ! Martix of D parameters for one dimesional foureir series
1571 do l=1,ntermd_1(i,j,k,iblock)
1572 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1573 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1574 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1575 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1576 ! write(iout,*) "whcodze" ,
1577 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1579 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1580 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1581 v2s(m,l,i,j,k,iblock),&
1582 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1583 ! Martix of D parameters for two dimesional fourier series
1584 do l=1,ntermd_2(i,j,k,iblock)
1586 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1587 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1588 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1589 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1598 write (iout,*) 'Constants for double torsionals'
1601 do j=-ntortyp+1,ntortyp-1
1602 do k=-ntortyp+1,ntortyp-1
1603 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1604 ' nsingle',ntermd_1(i,j,k,iblock),&
1605 ' ndouble',ntermd_2(i,j,k,iblock)
1607 write (iout,*) 'Single angles:'
1608 do l=1,ntermd_1(i,j,k,iblock)
1609 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1610 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1611 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1612 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1615 write (iout,*) 'Pairs of angles:'
1616 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1617 do l=1,ntermd_2(i,j,k,iblock)
1618 write (iout,'(i5,20f10.5)') &
1619 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1622 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1623 do l=1,ntermd_2(i,j,k,iblock)
1624 write (iout,'(i5,20f10.5)') &
1625 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1626 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1635 ! Read of Side-chain backbone correlation parameters
1636 ! Modified 11 May 2012 by Adasko
1639 read (isccor,*,end=119,err=119) nsccortyp
1641 !el from module energy-------------
1642 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1643 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1644 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1645 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1646 !-----------------------------------
1648 !el from module energy-------------
1649 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1651 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1653 isccortyp(i)=-isccortyp(-i)
1655 iscprol=isccortyp(20)
1656 ! write (iout,*) 'ntortyp',ntortyp
1658 !c maxinter is maximum interaction sites
1659 !el from module energy---------
1660 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1661 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1662 -nsccortyp:nsccortyp))
1663 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1664 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1665 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1666 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1667 !-----------------------------------
1671 read (isccor,*,end=119,err=119) &
1672 nterm_sccor(i,j),nlor_sccor(i,j)
1678 nterm_sccor(-i,j)=nterm_sccor(i,j)
1679 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1680 nterm_sccor(i,-j)=nterm_sccor(i,j)
1681 do k=1,nterm_sccor(i,j)
1682 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1684 if (j.eq.iscprol) then
1685 if (i.eq.isccortyp(10)) then
1686 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1687 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1689 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1690 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1691 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1692 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1693 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1694 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1695 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1696 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1699 if (i.eq.isccortyp(10)) then
1700 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1701 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1703 if (j.eq.isccortyp(10)) then
1704 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1705 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1707 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1708 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1709 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1710 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1711 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1712 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1716 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1717 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1718 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1719 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1722 do k=1,nlor_sccor(i,j)
1723 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1724 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1725 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1726 (1+vlor3sccor(k,i,j)**2)
1728 v0sccor(l,i,j)=v0ijsccor
1729 v0sccor(l,-i,j)=v0ijsccor1
1730 v0sccor(l,i,-j)=v0ijsccor2
1731 v0sccor(l,-i,-j)=v0ijsccor3
1737 !el from module energy-------------
1738 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1740 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1741 ! write (iout,*) 'ntortyp',ntortyp
1743 !c maxinter is maximum interaction sites
1744 !el from module energy---------
1745 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1746 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1747 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1748 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1749 !-----------------------------------
1753 read (isccor,*,end=119,err=119) &
1754 nterm_sccor(i,j),nlor_sccor(i,j)
1758 do k=1,nterm_sccor(i,j)
1759 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1761 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1764 do k=1,nlor_sccor(i,j)
1765 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1766 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1767 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1768 (1+vlor3sccor(k,i,j)**2)
1770 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1778 write (iout,'(/a/)') 'Torsional constants:'
1781 write (iout,*) 'ityp',i,' jtyp',j
1782 write (iout,*) 'Fourier constants'
1783 do k=1,nterm_sccor(i,j)
1784 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1786 write (iout,*) 'Lorenz constants'
1787 do k=1,nlor_sccor(i,j)
1788 write (iout,'(3(1pe15.5))') &
1789 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1796 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1797 ! interaction energy of the Gly, Ala, and Pro prototypes.
1802 write (iout,*) "Coefficients of the cumulants"
1804 read (ifourier,*) nloctyp
1805 !write(iout,*) "nloctyp",nloctyp
1806 !el from module energy-------
1807 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1808 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1809 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1810 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1811 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1812 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1813 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1814 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1818 !--------------------------------
1821 read (ifourier,*,end=115,err=115)
1822 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1824 write (iout,*) 'Type',i
1825 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1835 B1tilde(1,-i) =-b(3)
1837 ! b1tilde(1,i)=0.0d0
1838 ! b1tilde(2,i)=0.0d0
1863 Ctilde(1,2,-i)=-b(9)
1867 ! Ctilde(1,1,i)=0.0d0
1868 ! Ctilde(1,2,i)=0.0d0
1869 ! Ctilde(2,1,i)=0.0d0
1870 ! Ctilde(2,2,i)=0.0d0
1888 Dtilde(1,2,-i)=-b(8)
1892 ! Dtilde(1,1,i)=0.0d0
1893 ! Dtilde(1,2,i)=0.0d0
1894 ! Dtilde(2,1,i)=0.0d0
1895 ! Dtilde(2,2,i)=0.0d0
1896 EE(1,1,i)= b(10)+b(11)
1897 EE(2,2,i)=-b(10)+b(11)
1898 EE(2,1,i)= b(12)-b(13)
1899 EE(1,2,i)= b(12)+b(13)
1900 EE(1,1,-i)= b(10)+b(11)
1901 EE(2,2,-i)=-b(10)+b(11)
1902 EE(2,1,-i)=-b(12)+b(13)
1903 EE(1,2,-i)=-b(12)-b(13)
1909 ! ee(2,1,i)=ee(1,2,i)
1913 write (iout,*) 'Type',i
1915 write(iout,*) B1(1,i),B1(2,i)
1917 write(iout,*) B2(1,i),B2(2,i)
1920 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1924 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1928 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1933 ! Read electrostatic-interaction parameters
1938 write (iout,'(/a)') 'Electrostatic interaction constants:'
1939 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1940 'IT','JT','APP','BPP','AEL6','AEL3'
1942 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1943 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1944 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1945 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1950 app (i,j)=epp(i,j)*rri*rri
1951 bpp (i,j)=-2.0D0*epp(i,j)*rri
1952 ael6(i,j)=elpp6(i,j)*4.2D0**6
1953 ael3(i,j)=elpp3(i,j)*4.2D0**3
1955 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1961 ! Read side-chain interaction parameters.
1963 !el from module energy - COMMON.INTERACT-------
1964 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1965 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1966 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1967 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1968 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1979 !--------------------------------
1981 read (isidep,*,end=117,err=117) ipot,expon
1982 if (ipot.lt.1 .or. ipot.gt.5) then
1983 write (iout,'(2a)') 'Error while reading SC interaction',&
1984 'potential file - unknown potential type.'
1986 call MPI_Finalize(Ierror)
1992 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1993 ', exponents are ',expon,2*expon
1994 ! goto (10,20,30,30,40) ipot
1996 !----------------------- LJ potential ---------------------------------
1998 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1999 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2000 (sigma0(i),i=1,ntyp)
2002 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2003 write (iout,'(a/)') 'The epsilon array:'
2004 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2005 write (iout,'(/a)') 'One-body parameters:'
2006 write (iout,'(a,4x,a)') 'residue','sigma'
2007 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
2010 !----------------------- LJK potential --------------------------------
2012 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2013 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2014 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2016 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2017 write (iout,'(a/)') 'The epsilon array:'
2018 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2019 write (iout,'(/a)') 'One-body parameters:'
2020 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2021 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
2025 !---------------------- GB or BP potential -----------------------------
2029 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2031 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2032 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2033 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2034 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2035 ! For the GB potential convert sigma'**2 into chi'
2038 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2042 write (iout,'(/a/)') 'Parameters of the BP potential:'
2043 write (iout,'(a/)') 'The epsilon array:'
2044 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2045 write (iout,'(/a)') 'One-body parameters:'
2046 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2048 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
2049 chip(i),alp(i),i=1,ntyp)
2052 !--------------------- GBV potential -----------------------------------
2054 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2055 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2056 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2057 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2059 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2060 write (iout,'(a/)') 'The epsilon array:'
2061 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2062 write (iout,'(/a)') 'One-body parameters:'
2063 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2064 's||/s_|_^2',' chip ',' alph '
2065 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2066 sigii(i),chip(i),alp(i),i=1,ntyp)
2069 write(iout,*)"Wrong ipot"
2074 !-----------------------------------------------------------------------
2075 ! Calculate the "working" parameters of SC interactions.
2077 !el from module energy - COMMON.INTERACT-------
2078 allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2079 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2089 !--------------------------------
2098 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2099 sigma(j,i)=sigma(i,j)
2100 rs0(i,j)=dwa16*sigma(i,j)
2104 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2105 'Working parameters of the SC interactions:',&
2106 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2111 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2120 sigeps=dsign(1.0D0,epsij)
2122 aa(i,j)=epsij*rrij*rrij
2123 bb(i,j)=-sigeps*epsij*rrij
2127 sigt1sq=sigma0(i)**2
2128 sigt2sq=sigma0(j)**2
2131 ratsig1=sigt2sq/sigt1sq
2132 ratsig2=1.0D0/ratsig1
2133 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2134 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2135 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2139 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2140 sigmaii(i,j)=rsum_max
2141 sigmaii(j,i)=rsum_max
2143 ! sigmaii(i,j)=r0(i,j)
2144 ! sigmaii(j,i)=r0(i,j)
2146 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2147 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2148 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2149 augm(i,j)=epsij*r_augm**(2*expon)
2150 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2157 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2158 restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
2159 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2164 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2173 ! Define the SC-p interaction constants (hard-coded; old style)
2176 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2178 ! aad(i,1)=0.3D0*4.0D0**12
2179 ! Following line for constants currently implemented
2180 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2181 aad(i,1)=1.5D0*4.0D0**12
2182 ! aad(i,1)=0.17D0*5.6D0**12
2184 ! "Soft" SC-p repulsion
2186 ! Following line for constants currently implemented
2187 ! aad(i,1)=0.3D0*4.0D0**6
2188 ! "Hard" SC-p repulsion
2189 bad(i,1)=3.0D0*4.0D0**6
2190 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2199 ! 8/9/01 Read the SC-p interaction constants from file
2202 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2205 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2206 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2207 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2208 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2212 write (iout,*) "Parameters of SC-p interactions:"
2214 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2215 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2221 ! Define the constants of the disulfide bridge
2225 ! Old arbitrary potential - commented out.
2230 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2231 ! energy surface of diethyl disulfide.
2232 ! A. Liwo and U. Kozlowska, 11/24/03
2249 write (iout,'(/a)') "Disulfide bridge parameters:"
2250 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2251 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2252 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2253 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2257 111 write (iout,*) "Error reading bending energy parameters."
2259 112 write (iout,*) "Error reading rotamer energy parameters."
2261 113 write (iout,*) "Error reading torsional energy parameters."
2263 114 write (iout,*) "Error reading double torsional energy parameters."
2265 115 write (iout,*) &
2266 "Error reading cumulant (multibody energy) parameters."
2268 116 write (iout,*) "Error reading electrostatic energy parameters."
2270 117 write (iout,*) "Error reading side chain interaction parameters."
2272 118 write (iout,*) "Error reading SCp interaction parameters."
2274 119 write (iout,*) "Error reading SCCOR parameters"
2277 call MPI_Finalize(Ierror)
2281 end subroutine parmread
2283 !-----------------------------------------------------------------------------
2285 !-----------------------------------------------------------------------------
2286 subroutine printmat(ldim,m,n,iout,key,a)
2289 character(len=3),dimension(n) :: key
2290 real(kind=8),dimension(ldim,n) :: a
2292 integer :: i,j,k,m,iout,nlim
2296 write (iout,1000) (key(k),k=i,nlim)
2298 1000 format (/5x,8(6x,a3))
2299 1020 format (/80(1h-)/)
2301 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2304 1010 format (a3,2x,8(f9.4))
2306 end subroutine printmat
2307 !-----------------------------------------------------------------------------
2309 !-----------------------------------------------------------------------------
2311 ! Read the PDB file and convert the peptide geometry into virtual-chain
2314 use energy_data, only: itype
2318 use control, only: rescode
2319 ! implicit real*8 (a-h,o-z)
2320 ! include 'DIMENSIONS'
2321 ! include 'COMMON.LOCAL'
2322 ! include 'COMMON.VAR'
2323 ! include 'COMMON.CHAIN'
2324 ! include 'COMMON.INTERACT'
2325 ! include 'COMMON.IOUNITS'
2326 ! include 'COMMON.GEO'
2327 ! include 'COMMON.NAMES'
2328 ! include 'COMMON.CONTROL'
2329 ! include 'COMMON.DISTFIT'
2330 ! include 'COMMON.SETUP'
2331 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2333 logical :: lprn=.true.,fail
2334 real(kind=8),dimension(3) :: e1,e2,e3
2335 real(kind=8) :: dcj,efree_temp
2336 character(len=3) :: seq,res
2337 character(len=5) :: atom
2338 character(len=80) :: card
2339 real(kind=8),dimension(3,20) :: sccor
2340 integer :: kkk,lll,icha,kupa !rescode,
2343 integer,dimension(2,maxres/3) :: hfrag_alloc
2344 integer,dimension(4,maxres/3) :: bfrag_alloc
2345 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2351 ! write (2,*) "UNRES_PDB",unres_pdb
2359 !-----------------------------
2360 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2361 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2364 read (ipdbin,'(a80)',end=10) card
2365 ! write (iout,'(a)') card
2366 if (card(:5).eq.'HELIX') then
2369 read(card(22:25),*) hfrag(1,nhfrag)
2370 read(card(34:37),*) hfrag(2,nhfrag)
2372 if (card(:5).eq.'SHEET') then
2375 read(card(24:26),*) bfrag(1,nbfrag)
2376 read(card(35:37),*) bfrag(2,nbfrag)
2377 !rc----------------------------------------
2378 !rc to be corrected !!!
2379 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2380 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2381 !rc----------------------------------------
2383 if (card(:3).eq.'END') then
2385 else if (card(:3).eq.'TER') then
2389 itype(ires_old)=ntyp1
2391 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2394 dc(j,ires)=sccor(j,iii)
2397 call sccenter(ires,iii,sccor)
2402 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2403 ! Fish out the ATOM cards.
2404 if (index(card(1:4),'ATOM').gt.0) then
2405 read (card(12:16),*) atom
2406 ! write (iout,*) "! ",atom," !",ires
2407 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2408 read (card(23:26),*) ires
2409 read (card(18:20),'(a3)') res
2410 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2411 ! & " ires_old",ires_old
2412 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2413 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2414 if (ires-ishift+ishift1.ne.ires_old) then
2415 ! Calculate the CM of the preceding residue.
2416 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2418 ! write (iout,*) "Calculating sidechain center iii",iii
2421 dc(j,ires+nres)=sccor(j,iii)
2424 call sccenter(ires_old,iii,sccor)
2428 ! Start new residue.
2429 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2432 else if (ibeg.eq.1) then
2433 write (iout,*) "BEG ires",ires
2435 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2439 ires=ires-ishift+ishift1
2441 ! write (iout,*) "ishift",ishift," ires",ires,&
2442 ! " ires_old",ires_old
2444 else if (ibeg.eq.2) then
2446 ishift=-ires_old+ires-1 !!!!!
2447 ishift1=ishift1-1 !!!!!
2448 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2449 ires=ires-ishift+ishift1
2453 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2454 ires=ires-ishift+ishift1
2457 if (res.eq.'ACE' .or. res.eq.'NHE') then
2460 itype(ires)=rescode(ires,res,0)
2463 ires=ires-ishift+ishift1
2465 ! write (iout,*) "ires_old",ires_old," ires",ires
2466 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2469 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2470 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2471 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2472 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2473 ! write (iout,*) "backbone ",atom
2475 write (iout,'(2i3,2x,a,3f8.3)') &
2476 ires,itype(ires),res,(c(j,ires),j=1,3)
2480 sccor(j,iii)=c(j,ires)
2482 ! write (*,*) card(23:27),ires,itype(ires)
2483 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2484 atom.ne.'N' .and. atom.ne.'C' .and. &
2485 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2486 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2487 ! write (iout,*) "sidechain ",atom
2489 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2493 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2494 if (ires.eq.0) return
2495 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2499 ! write (iout,*) i,itype(i)
2500 if (itype(i).eq.ntyp1) then
2501 ! write (iout,*) "dummy",i,itype(i)
2503 c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2504 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2509 ! Calculate the CM of the last side chain.
2513 dc(j,ires)=sccor(j,iii)
2516 call sccenter(ires,iii,sccor)
2522 if (itype(nres).ne.10) then
2526 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2527 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2534 c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
2538 dcj=c(j,nres-2)-c(j,nres-3)
2539 c(j,nres)=c(j,nres-1)+dcj
2540 c(j,2*nres)=c(j,nres)
2544 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2546 if (nres.ne.nres0) then
2547 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2549 stop "Error nres value in WHAM input"
2552 !---------------------------------
2553 !el reallocate tables
2556 ! hfrag_alloc(j,i)=hfrag(j,i)
2559 ! bfrag_alloc(j,i)=bfrag(j,i)
2565 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2566 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2567 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2568 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2572 ! hfrag(j,i)=hfrag_alloc(j,i)
2577 ! bfrag(j,i)=bfrag_alloc(j,i)
2580 !el end reallocate tables
2581 !---------------------------------
2589 c(j,2*nres)=c(j,nres)
2591 if (itype(1).eq.ntyp1) then
2595 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2596 call refsys(2,3,4,e1,e2,e3,fail)
2603 c(j,1)=c(j,2)-3.8d0*e2(j)
2613 ! Copy the coordinates to reference coordinates
2619 ! Calculate internal coordinates.
2621 write (iout,'(/a)') &
2622 "Cartesian coordinates of the reference structure"
2623 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2624 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2626 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2627 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2628 (c(j,ires+nres),j=1,3)
2631 ! znamy już nres wiec mozna alokowac tablice
2632 ! Calculate internal coordinates.
2633 if(me.eq.king.or..not.out1file)then
2634 write (iout,'(a)') &
2635 "Backbone and SC coordinates as read from the PDB"
2637 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2638 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2639 (c(j,nres+ires),j=1,3)
2643 if(.not.allocated(vbld)) then
2644 allocate(vbld(2*nres))
2649 if(.not.allocated(vbld_inv)) then
2650 allocate(vbld_inv(2*nres))
2656 if(.not.allocated(theta)) then
2657 allocate(theta(nres+2))
2658 ! allocate(phi(nres+2))
2659 ! allocate(alph(nres+2))
2660 ! allocate(omeg(nres+2))
2668 ! allocate(costtab(nres))
2669 ! allocate(sinttab(nres))
2670 ! allocate(cost2tab(nres))
2671 ! allocate(sint2tab(nres))
2672 ! allocate(xxref(nres))
2673 ! allocate(yyref(nres))
2674 ! allocate(zzref(nres)) !(maxres)
2686 if(.not.allocated(phi)) allocate(phi(nres+2))
2687 if(.not.allocated(alph)) allocate(alph(nres+2))
2688 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2689 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2690 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2691 if(.not.allocated(costtab)) allocate(costtab(nres))
2692 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2693 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2694 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2695 if(.not.allocated(xxref)) allocate(xxref(nres))
2696 if(.not.allocated(yyref)) allocate(yyref(nres))
2697 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2698 if(.not.allocated(dc_norm)) then
2699 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2700 allocate(dc_norm(3,0:2*nres+2))
2708 call int_from_cart(.true.,.false.)
2709 call sc_loc_geom(.true.)
2710 ! call sc_loc_geom(.false.)
2711 ! wczesbiej bylo false
2713 thetaref(i)=theta(i)
2723 dc(j,i)=c(j,i+1)-c(j,i)
2724 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2729 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2730 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2732 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2736 ! Copy the coordinates to reference coordinates
2737 ! Splits to single chain if occurs
2741 ! cref(j,i,cou)=c(j,i)
2745 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2746 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2747 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2748 !-----------------------------
2754 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2756 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2759 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2764 cref(j,i,cou)=c(j,i)
2765 cref(j,i+nres,cou)=c(j,i+nres)
2767 chain_rep(j,lll,kkk)=c(j,i)
2768 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2772 write (iout,*) chain_length
2773 if (chain_length.eq.0) chain_length=nres
2775 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2776 chain_rep(j,chain_length+nres,symetr) &
2777 =chain_rep(j,chain_length+nres,1)
2780 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2782 ! do kkk=1,chain_length
2783 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2787 ! makes copy of chains
2788 write (iout,*) "symetr", symetr
2790 if (symetr.gt.1) then
2797 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2803 ! write (iout,*) i,icha
2804 do lll=1,chain_length
2806 if (cou.le.nres) then
2808 kupa=mod(lll,chain_length)
2809 iprzes=(kkk-1)*chain_length+lll
2810 if (kupa.eq.0) kupa=chain_length
2811 ! write (iout,*) "kupa", kupa
2812 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2813 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2820 !-koniec robienia kopii
2823 write (iout,*) "nowa struktura", nperm
2825 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2827 cref(3,i,kkk),cref(1,nres+i,kkk),&
2828 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2830 100 format (//' alpha-carbon coordinates ',&
2831 ' centroid coordinates'/ &
2832 ' ', 6X,'X',11X,'Y',11X,'Z', &
2833 10X,'X',11X,'Y',11X,'Z')
2834 110 format (a,'(',i3,')',6f12.5)
2840 bfrag(i,j)=bfrag(i,j)-ishift
2846 hfrag(i,j)=hfrag(i,j)-ishift
2852 end subroutine readpdb
2853 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2854 !-----------------------------------------------------------------------------
2856 !-----------------------------------------------------------------------------
2857 subroutine read_control
2871 use random, only: random_init
2872 ! implicit real*8 (a-h,o-z)
2873 ! include 'DIMENSIONS'
2875 use prng, only:prng_restart
2877 logical :: OKRandom!, prng_restart
2880 ! include 'COMMON.IOUNITS'
2881 ! include 'COMMON.TIME1'
2882 ! include 'COMMON.THREAD'
2883 ! include 'COMMON.SBRIDGE'
2884 ! include 'COMMON.CONTROL'
2885 ! include 'COMMON.MCM'
2886 ! include 'COMMON.MAP'
2887 ! include 'COMMON.HEADER'
2888 ! include 'COMMON.CSA'
2889 ! include 'COMMON.CHAIN'
2890 ! include 'COMMON.MUCA'
2891 ! include 'COMMON.MD'
2892 ! include 'COMMON.FFIELD'
2893 ! include 'COMMON.INTERACT'
2894 ! include 'COMMON.SETUP'
2895 !el integer :: KDIAG,ICORFL,IXDR
2896 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2897 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2898 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2899 ! character(len=80) :: ucase
2900 character(len=640) :: controlcard
2902 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2908 read (INP,'(a)') titel
2909 call card_concat(controlcard,.true.)
2910 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2911 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2912 call reada(controlcard,'SEED',seed,0.0D0)
2913 call random_init(seed)
2914 ! Set up the time limit (caution! The time must be input in minutes!)
2915 read_cart=index(controlcard,'READ_CART').gt.0
2916 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2917 call readi(controlcard,'SYM',symetr,1)
2918 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2919 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2920 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2921 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2922 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2923 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2924 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2925 call reada(controlcard,'DRMS',drms,0.1D0)
2926 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2927 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2928 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2929 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2930 write (iout,'(a,f10.1)')'DRMS = ',drms
2931 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2932 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2934 call readi(controlcard,'NZ_START',nz_start,0)
2935 call readi(controlcard,'NZ_END',nz_end,0)
2936 call readi(controlcard,'IZ_SC',iz_sc,0)
2937 timlim=60.0D0*timlim
2938 safety = 60.0d0*safety
2941 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2942 minim=(index(controlcard,'MINIMIZE').gt.0)
2943 dccart=(index(controlcard,'CART').gt.0)
2944 overlapsc=(index(controlcard,'OVERLAP').gt.0)
2945 overlapsc=.not.overlapsc
2946 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2947 searchsc=.not.searchsc
2948 sideadd=(index(controlcard,'SIDEADD').gt.0)
2949 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2950 outpdb=(index(controlcard,'PDBOUT').gt.0)
2951 outmol2=(index(controlcard,'MOL2OUT').gt.0)
2952 pdbref=(index(controlcard,'PDBREF').gt.0)
2953 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2954 indpdb=index(controlcard,'PDBSTART')
2955 extconf=(index(controlcard,'EXTCONF').gt.0)
2956 call readi(controlcard,'IPRINT',iprint,0)
2957 call readi(controlcard,'MAXGEN',maxgen,10000)
2958 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2959 call readi(controlcard,"KDIAG",kdiag,0)
2960 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2961 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2962 write (iout,*) "RESCALE_MODE",rescale_mode
2963 split_ene=index(controlcard,'SPLIT_ENE').gt.0
2964 if (index(controlcard,'REGULAR').gt.0.0D0) then
2965 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2969 if (index(controlcard,'CHECKGRAD').gt.0) then
2971 if (index(controlcard,'CART').gt.0) then
2973 elseif (index(controlcard,'CARINT').gt.0) then
2978 elseif (index(controlcard,'THREAD').gt.0) then
2980 call readi(controlcard,'THREAD',nthread,0)
2981 if (nthread.gt.0) then
2982 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2985 write (iout,'(a)')'A number has to follow the THREAD keyword.'
2986 stop 'Error termination in Read_Control.'
2988 else if (index(controlcard,'MCMA').gt.0) then
2990 else if (index(controlcard,'MCEE').gt.0) then
2992 else if (index(controlcard,'MULTCONF').gt.0) then
2994 else if (index(controlcard,'MAP').gt.0) then
2996 call readi(controlcard,'MAP',nmap,0)
2997 else if (index(controlcard,'CSA').gt.0) then
2999 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3001 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3004 !fcm else if (index(controlcard,'MCMF').gt.0) then
3006 else if (index(controlcard,'SOFTREG').gt.0) then
3008 else if (index(controlcard,'CHECK_BOND').gt.0) then
3010 else if (index(controlcard,'TEST').gt.0) then
3012 else if (index(controlcard,'MD').gt.0) then
3014 else if (index(controlcard,'RE ').gt.0) then
3018 lmuca=index(controlcard,'MUCA').gt.0
3019 call readi(controlcard,'MUCADYN',mucadyn,0)
3020 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3021 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3023 write (iout,*) 'MUCADYN=',mucadyn
3024 write (iout,*) 'MUCASMOOTH=',muca_smooth
3027 iscode=index(controlcard,'ONE_LETTER')
3028 indphi=index(controlcard,'PHI')
3029 indback=index(controlcard,'BACK')
3030 iranconf=index(controlcard,'RAND_CONF')
3031 i2ndstr=index(controlcard,'USE_SEC_PRED')
3032 gradout=index(controlcard,'GRADOUT').gt.0
3033 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3034 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3035 if (me.eq.king .or. .not.out1file ) &
3036 write (iout,*) "DISTCHAINMAX",distchainmax
3038 if(me.eq.king.or..not.out1file) &
3039 write (iout,'(2a)') diagmeth(kdiag),&
3040 ' routine used to diagonalize matrices.'
3042 end subroutine read_control
3043 !-----------------------------------------------------------------------------
3044 subroutine read_REMDpar
3046 ! Read REMD settings
3053 use control_data, only:out1file
3054 ! implicit real*8 (a-h,o-z)
3055 ! include 'DIMENSIONS'
3056 ! include 'COMMON.IOUNITS'
3057 ! include 'COMMON.TIME1'
3058 ! include 'COMMON.MD'
3061 !el include 'COMMON.LANGEVIN'
3063 !el include 'COMMON.LANGEVIN.lang0'
3065 ! include 'COMMON.INTERACT'
3066 ! include 'COMMON.NAMES'
3067 ! include 'COMMON.GEO'
3068 ! include 'COMMON.REMD'
3069 ! include 'COMMON.CONTROL'
3070 ! include 'COMMON.SETUP'
3071 ! character(len=80) :: ucase
3072 character(len=320) :: controlcard
3073 character(len=3200) :: controlcard1
3074 integer :: iremd_m_total
3077 ! real(kind=8) :: var,ene
3079 if(me.eq.king.or..not.out1file) &
3080 write (iout,*) "REMD setup"
3082 call card_concat(controlcard,.true.)
3083 call readi(controlcard,"NREP",nrep,3)
3084 call readi(controlcard,"NSTEX",nstex,1000)
3085 call reada(controlcard,"RETMIN",retmin,10.0d0)
3086 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3087 mremdsync=(index(controlcard,'SYNC').gt.0)
3088 call readi(controlcard,"NSYN",i_sync_step,100)
3089 restart1file=(index(controlcard,'REST1FILE').gt.0)
3090 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3091 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3092 if(max_cache_traj_use.gt.max_cache_traj) &
3093 max_cache_traj_use=max_cache_traj
3094 if(me.eq.king.or..not.out1file) then
3095 !d if (traj1file) then
3096 !rc caching is in testing - NTWX is not ignored
3097 !d write (iout,*) "NTWX value is ignored"
3098 !d write (iout,*) " trajectory is stored to one file by master"
3099 !d write (iout,*) " before exchange at NSTEX intervals"
3101 write (iout,*) "NREP= ",nrep
3102 write (iout,*) "NSTEX= ",nstex
3103 write (iout,*) "SYNC= ",mremdsync
3104 write (iout,*) "NSYN= ",i_sync_step
3105 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3108 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3109 if (index(controlcard,'TLIST').gt.0) then
3111 call card_concat(controlcard1,.true.)
3112 read(controlcard1,*) (remd_t(i),i=1,nrep)
3113 if(me.eq.king.or..not.out1file) &
3114 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3117 if (index(controlcard,'MLIST').gt.0) then
3119 call card_concat(controlcard1,.true.)
3120 read(controlcard1,*) (remd_m(i),i=1,nrep)
3121 if(me.eq.king.or..not.out1file) then
3122 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3125 iremd_m_total=iremd_m_total+remd_m(i)
3127 write (iout,*) 'Total number of replicas ',iremd_m_total
3130 if(me.eq.king.or..not.out1file) &
3131 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3133 end subroutine read_REMDpar
3134 !-----------------------------------------------------------------------------
3135 subroutine read_MDpar
3139 use control_data, only: r_cut,rlamb,out1file
3141 use geometry_data, only: pi
3143 ! implicit real*8 (a-h,o-z)
3144 ! include 'DIMENSIONS'
3145 ! include 'COMMON.IOUNITS'
3146 ! include 'COMMON.TIME1'
3147 ! include 'COMMON.MD'
3150 !el include 'COMMON.LANGEVIN'
3152 !el include 'COMMON.LANGEVIN.lang0'
3154 ! include 'COMMON.INTERACT'
3155 ! include 'COMMON.NAMES'
3156 ! include 'COMMON.GEO'
3157 ! include 'COMMON.SETUP'
3158 ! include 'COMMON.CONTROL'
3159 ! include 'COMMON.SPLITELE'
3160 ! character(len=80) :: ucase
3161 character(len=320) :: controlcard
3166 call card_concat(controlcard,.true.)
3167 call readi(controlcard,"NSTEP",n_timestep,1000000)
3168 call readi(controlcard,"NTWE",ntwe,100)
3169 call readi(controlcard,"NTWX",ntwx,1000)
3170 call reada(controlcard,"DT",d_time,1.0d-1)
3171 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3172 call reada(controlcard,"DAMAX",damax,1.0d1)
3173 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3174 call readi(controlcard,"LANG",lang,0)
3175 RESPA = index(controlcard,"RESPA") .gt. 0
3176 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3177 ntime_split0=ntime_split
3178 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3179 ntime_split0=ntime_split
3180 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3181 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3182 rest = index(controlcard,"REST").gt.0
3183 tbf = index(controlcard,"TBF").gt.0
3184 usampl = index(controlcard,"USAMPL").gt.0
3185 mdpdb = index(controlcard,"MDPDB").gt.0
3186 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3187 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3188 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3189 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3190 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3191 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3192 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3193 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3194 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3195 large = index(controlcard,"LARGE").gt.0
3196 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3197 rattle = index(controlcard,"RATTLE").gt.0
3198 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3204 if(me.eq.king.or..not.out1file) then
3206 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3208 write (iout,'(a)') "The units are:"
3209 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3210 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3211 " acceleration: angstrom/(48.9 fs)**2"
3212 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3214 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3215 write (iout,'(a60,f10.5,a)') &
3216 "Initial time step of numerical integration:",d_time,&
3218 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3220 write (iout,'(2a,i4,a)') &
3221 "A-MTS algorithm used; initial time step for fast-varying",&
3222 " short-range forces split into",ntime_split," steps."
3223 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3224 r_cut," lambda",rlamb
3226 write (iout,'(2a,f10.5)') &
3227 "Maximum acceleration threshold to reduce the time step",&
3228 "/increase split number:",damax
3229 write (iout,'(2a,f10.5)') &
3230 "Maximum predicted energy drift to reduce the timestep",&
3231 "/increase split number:",edriftmax
3232 write (iout,'(a60,f10.5)') &
3233 "Maximum velocity threshold to reduce velocities:",dvmax
3234 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3235 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3236 if (rattle) write (iout,'(a60)') &
3237 "Rattle algorithm used to constrain the virtual bonds"
3241 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3242 call reada(controlcard,"RWAT",rwat,1.4d0)
3243 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3244 surfarea=index(controlcard,"SURFAREA").gt.0
3245 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3246 if(me.eq.king.or..not.out1file)then
3247 write (iout,'(/a,$)') "Langevin dynamics calculation"
3249 write (iout,'(a/)') &
3250 " with direct integration of Langevin equations"
3251 else if (lang.eq.2) then
3252 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3253 else if (lang.eq.3) then
3254 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3255 else if (lang.eq.4) then
3256 write (iout,'(a/)') " in overdamped mode"
3258 write (iout,'(//a,i5)') &
3259 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3262 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3263 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3264 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3265 write (iout,'(a60,f10.5)') &
3266 "Scaling factor of the friction forces:",scal_fric
3267 if (surfarea) write (iout,'(2a,i10,a)') &
3268 "Friction coefficients will be scaled by solvent-accessible",&
3269 " surface area every",reset_fricmat," steps."
3271 ! Calculate friction coefficients and bounds of stochastic forces
3272 eta=6*pi*cPoise*etawat
3273 if(me.eq.king.or..not.out1file) &
3274 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3276 gamp=scal_fric*(pstok+rwat)*eta
3277 stdfp=dsqrt(2*Rb*t_bath/d_time)
3278 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3280 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3281 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3283 if(me.eq.king.or..not.out1file)then
3284 write (iout,'(/2a/)') &
3285 "Radii of site types and friction coefficients and std's of",&
3286 " stochastic forces of fully exposed sites"
3287 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3289 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3290 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3294 if(me.eq.king.or..not.out1file)then
3295 write (iout,'(a)') "Berendsen bath calculation"
3296 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3297 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3299 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3300 count_reset_moment," steps"
3302 write (iout,'(a,i10,a)') &
3303 "Velocities will be reset at random every",count_reset_vel,&
3307 if(me.eq.king.or..not.out1file) &
3308 write (iout,'(a31)') "Microcanonical mode calculation"
3310 if(me.eq.king.or..not.out1file)then
3311 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3313 write(iout,*) "MD running with constraints."
3314 write(iout,*) "Equilibration time ", eq_time, " mtus."
3315 write(iout,*) "Constraining ", nfrag," fragments."
3316 write(iout,*) "Length of each fragment, weight and q0:"
3318 write (iout,*) "Set of restraints #",iset
3320 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3321 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3323 write(iout,*) "constraints between ", npair, "fragments."
3324 write(iout,*) "constraint pairs, weights and q0:"
3326 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3327 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3329 write(iout,*) "angle constraints within ", nfrag_back,&
3330 "backbone fragments."
3331 write(iout,*) "fragment, weights:"
3333 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3334 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3335 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3338 iset=mod(kolor,nset)+1
3341 if(me.eq.king.or..not.out1file) &
3342 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3344 end subroutine read_MDpar
3345 !-----------------------------------------------------------------------------
3349 ! implicit real*8 (a-h,o-z)
3350 ! include 'DIMENSIONS'
3351 ! include 'COMMON.MAP'
3352 ! include 'COMMON.IOUNITS'
3353 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3354 character(len=80) :: mapcard !,ucase
3357 ! real(kind=8) :: var,ene
3360 read (inp,'(a)') mapcard
3361 mapcard=ucase(mapcard)
3362 if (index(mapcard,'PHI').gt.0) then
3364 else if (index(mapcard,'THE').gt.0) then
3366 else if (index(mapcard,'ALP').gt.0) then
3368 else if (index(mapcard,'OME').gt.0) then
3371 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3372 stop 'Error - illegal variable spec in MAP card.'
3374 call readi (mapcard,'RES1',res1(imap),0)
3375 call readi (mapcard,'RES2',res2(imap),0)
3376 if (res1(imap).eq.0) then
3377 res1(imap)=res2(imap)
3378 else if (res2(imap).eq.0) then
3379 res2(imap)=res1(imap)
3381 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3382 write (iout,'(a)') &
3383 'Error - illegal definition of variable group in MAP.'
3384 stop 'Error - illegal definition of variable group in MAP.'
3386 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3387 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3388 call readi(mapcard,'NSTEP',nstep(imap),0)
3389 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3390 write (iout,'(a)') &
3391 'Illegal boundary and/or step size specification in MAP.'
3392 stop 'Illegal boundary and/or step size specification in MAP.'
3396 end subroutine map_read
3397 !-----------------------------------------------------------------------------
3400 use control_data, only: vdisulf
3402 ! implicit real*8 (a-h,o-z)
3403 ! include 'DIMENSIONS'
3404 ! include 'COMMON.IOUNITS'
3405 ! include 'COMMON.GEO'
3406 ! include 'COMMON.CSA'
3407 ! include 'COMMON.BANK'
3408 ! include 'COMMON.CONTROL'
3409 ! character(len=80) :: ucase
3410 character(len=620) :: mcmcard
3412 ! integer :: ntf,ik,iw_pdb
3413 ! real(kind=8) :: var,ene
3415 call card_concat(mcmcard,.true.)
3417 call readi(mcmcard,'NCONF',nconf,50)
3418 call readi(mcmcard,'NADD',nadd,0)
3419 call readi(mcmcard,'JSTART',jstart,1)
3420 call readi(mcmcard,'JEND',jend,1)
3421 call readi(mcmcard,'NSTMAX',nstmax,500000)
3422 call readi(mcmcard,'N0',n0,1)
3423 call readi(mcmcard,'N1',n1,6)
3424 call readi(mcmcard,'N2',n2,4)
3425 call readi(mcmcard,'N3',n3,0)
3426 call readi(mcmcard,'N4',n4,0)
3427 call readi(mcmcard,'N5',n5,0)
3428 call readi(mcmcard,'N6',n6,10)
3429 call readi(mcmcard,'N7',n7,0)
3430 call readi(mcmcard,'N8',n8,0)
3431 call readi(mcmcard,'N9',n9,0)
3432 call readi(mcmcard,'N14',n14,0)
3433 call readi(mcmcard,'N15',n15,0)
3434 call readi(mcmcard,'N16',n16,0)
3435 call readi(mcmcard,'N17',n17,0)
3436 call readi(mcmcard,'N18',n18,0)
3438 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3440 call readi(mcmcard,'NDIFF',ndiff,2)
3441 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3442 call readi(mcmcard,'IS1',is1,1)
3443 call readi(mcmcard,'IS2',is2,8)
3444 call readi(mcmcard,'NRAN0',nran0,4)
3445 call readi(mcmcard,'NRAN1',nran1,2)
3446 call readi(mcmcard,'IRR',irr,1)
3447 call readi(mcmcard,'NSEED',nseed,20)
3448 call readi(mcmcard,'NTOTAL',ntotal,10000)
3449 call reada(mcmcard,'CUT1',cut1,2.0d0)
3450 call reada(mcmcard,'CUT2',cut2,5.0d0)
3451 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3452 call readi(mcmcard,'ICMAX',icmax,3)
3453 call readi(mcmcard,'IRESTART',irestart,0)
3454 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3457 call reada(mcmcard,'DELE',dele,20.0d0)
3458 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3459 call readi(mcmcard,'IREF',iref,0)
3460 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3461 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3462 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3463 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3464 write (iout,*) "NCONF_IN",nconf_in
3466 end subroutine csaread
3467 !-----------------------------------------------------------------------------
3471 use control_data, only: MaxMoveType
3474 ! implicit real*8 (a-h,o-z)
3475 ! include 'DIMENSIONS'
3476 ! include 'COMMON.MCM'
3477 ! include 'COMMON.MCE'
3478 ! include 'COMMON.IOUNITS'
3479 ! character(len=80) :: ucase
3480 character(len=320) :: mcmcard
3483 ! real(kind=8) :: var,ene
3485 call card_concat(mcmcard,.true.)
3486 call readi(mcmcard,'MAXACC',maxacc,100)
3487 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3488 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3489 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3490 call readi(mcmcard,'MAXREPM',maxrepm,200)
3491 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3492 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3493 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3494 call reada(mcmcard,'E_UP',e_up,5.0D0)
3495 call reada(mcmcard,'DELTE',delte,0.1D0)
3496 call readi(mcmcard,'NSWEEP',nsweep,5)
3497 call readi(mcmcard,'NSTEPH',nsteph,0)
3498 call readi(mcmcard,'NSTEPC',nstepc,0)
3499 call reada(mcmcard,'TMIN',tmin,298.0D0)
3500 call reada(mcmcard,'TMAX',tmax,298.0D0)
3501 call readi(mcmcard,'NWINDOW',nwindow,0)
3502 call readi(mcmcard,'PRINT_MC',print_mc,0)
3503 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3504 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3505 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3506 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3507 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3508 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3509 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3510 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3511 if (nwindow.gt.0) then
3512 allocate(winstart(nwindow)) !!el (maxres)
3513 allocate(winend(nwindow)) !!el
3514 allocate(winlen(nwindow)) !!el
3515 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3517 winlen(i)=winend(i)-winstart(i)+1
3520 if (tmax.lt.tmin) tmax=tmin
3521 if (tmax.eq.tmin) then
3525 if (nstepc.gt.0 .and. nsteph.gt.0) then
3526 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3527 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3529 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3530 ! Probabilities of different move types
3531 sumpro_type(0)=0.0D0
3532 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3533 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3534 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3535 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3536 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3537 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3538 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3540 print *,'i',i,' sumprotype',sumpro_type(i)
3541 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3542 print *,'i',i,' sumprotype',sumpro_type(i)
3545 end subroutine mcmread
3546 !-----------------------------------------------------------------------------
3547 subroutine read_minim
3550 ! implicit real*8 (a-h,o-z)
3551 ! include 'DIMENSIONS'
3552 ! include 'COMMON.MINIM'
3553 ! include 'COMMON.IOUNITS'
3554 ! character(len=80) :: ucase
3555 character(len=320) :: minimcard
3557 ! integer :: ntf,ik,iw_pdb
3558 ! real(kind=8) :: var,ene
3560 call card_concat(minimcard,.true.)
3561 call readi(minimcard,'MAXMIN',maxmin,2000)
3562 call readi(minimcard,'MAXFUN',maxfun,5000)
3563 call readi(minimcard,'MINMIN',minmin,maxmin)
3564 call readi(minimcard,'MINFUN',minfun,maxmin)
3565 call reada(minimcard,'TOLF',tolf,1.0D-2)
3566 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3567 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3568 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3569 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3570 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3571 'Options in energy minimization:'
3572 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3573 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3574 'MinMin:',MinMin,' MinFun:',MinFun,&
3575 ' TolF:',TolF,' RTolF:',RTolF
3577 end subroutine read_minim
3578 !-----------------------------------------------------------------------------
3579 subroutine openunits
3581 use energy_data, only: usampl
3584 use control_data, only:out1file
3585 use control, only: getenv_loc
3586 ! implicit real*8 (a-h,o-z)
3587 ! include 'DIMENSIONS'
3590 character(len=16) :: form,nodename
3591 integer :: nodelen,ierror,npos
3593 ! include 'COMMON.SETUP'
3594 ! include 'COMMON.IOUNITS'
3595 ! include 'COMMON.MD'
3596 ! include 'COMMON.CONTROL'
3597 integer :: lenpre,lenpot,lentmp !,ilen
3599 character(len=3) :: out1file_text !,ucase
3600 character(len=3) :: ll
3603 ! integer :: ntf,ik,iw_pdb
3604 ! real(kind=8) :: var,ene
3606 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3607 call getenv_loc("PREFIX",prefix)
3609 call getenv_loc("POT",pot)
3610 call getenv_loc("DIRTMP",tmpdir)
3611 call getenv_loc("CURDIR",curdir)
3612 call getenv_loc("OUT1FILE",out1file_text)
3613 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3614 out1file_text=ucase(out1file_text)
3615 if (out1file_text(1:1).eq."Y") then
3618 out1file=fg_rank.gt.0
3623 if (lentmp.gt.0) then
3624 write (*,'(80(1h!))')
3625 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3626 write (*,'(80(1h!))')
3627 write (*,*)"All output files will be on node /tmp directory."
3629 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3630 if (me.eq.king) then
3631 write (*,*) "The master node is ",nodename
3632 else if (fg_rank.eq.0) then
3633 write (*,*) "I am the CG slave node ",nodename
3635 write (*,*) "I am the FG slave node ",nodename
3638 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3639 lenpre = lentmp+lenpre+1
3641 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3642 ! Get the names and open the input files
3643 #if defined(WINIFL) || defined(WINPGI)
3644 open(1,file=pref_orig(:ilen(pref_orig))// &
3645 '.inp',status='old',readonly,shared)
3646 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3647 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3648 ! Get parameter filenames and open the parameter files.
3649 call getenv_loc('BONDPAR',bondname)
3650 open (ibond,file=bondname,status='old',readonly,shared)
3651 call getenv_loc('THETPAR',thetname)
3652 open (ithep,file=thetname,status='old',readonly,shared)
3653 call getenv_loc('ROTPAR',rotname)
3654 open (irotam,file=rotname,status='old',readonly,shared)
3655 call getenv_loc('TORPAR',torname)
3656 open (itorp,file=torname,status='old',readonly,shared)
3657 call getenv_loc('TORDPAR',tordname)
3658 open (itordp,file=tordname,status='old',readonly,shared)
3659 call getenv_loc('FOURIER',fouriername)
3660 open (ifourier,file=fouriername,status='old',readonly,shared)
3661 call getenv_loc('ELEPAR',elename)
3662 open (ielep,file=elename,status='old',readonly,shared)
3663 call getenv_loc('SIDEPAR',sidename)
3664 open (isidep,file=sidename,status='old',readonly,shared)
3665 #elif (defined CRAY) || (defined AIX)
3666 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3668 ! print *,"Processor",myrank," opened file 1"
3669 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3670 ! print *,"Processor",myrank," opened file 9"
3671 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3672 ! Get parameter filenames and open the parameter files.
3673 call getenv_loc('BONDPAR',bondname)
3674 open (ibond,file=bondname,status='old',action='read')
3675 ! print *,"Processor",myrank," opened file IBOND"
3676 call getenv_loc('THETPAR',thetname)
3677 open (ithep,file=thetname,status='old',action='read')
3678 ! print *,"Processor",myrank," opened file ITHEP"
3679 call getenv_loc('ROTPAR',rotname)
3680 open (irotam,file=rotname,status='old',action='read')
3681 ! print *,"Processor",myrank," opened file IROTAM"
3682 call getenv_loc('TORPAR',torname)
3683 open (itorp,file=torname,status='old',action='read')
3684 ! print *,"Processor",myrank," opened file ITORP"
3685 call getenv_loc('TORDPAR',tordname)
3686 open (itordp,file=tordname,status='old',action='read')
3687 ! print *,"Processor",myrank," opened file ITORDP"
3688 call getenv_loc('SCCORPAR',sccorname)
3689 open (isccor,file=sccorname,status='old',action='read')
3690 ! print *,"Processor",myrank," opened file ISCCOR"
3691 call getenv_loc('FOURIER',fouriername)
3692 open (ifourier,file=fouriername,status='old',action='read')
3693 ! print *,"Processor",myrank," opened file IFOURIER"
3694 call getenv_loc('ELEPAR',elename)
3695 open (ielep,file=elename,status='old',action='read')
3696 ! print *,"Processor",myrank," opened file IELEP"
3697 call getenv_loc('SIDEPAR',sidename)
3698 open (isidep,file=sidename,status='old',action='read')
3699 ! print *,"Processor",myrank," opened file ISIDEP"
3700 ! print *,"Processor",myrank," opened parameter files"
3702 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3703 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3704 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3705 ! Get parameter filenames and open the parameter files.
3706 call getenv_loc('BONDPAR',bondname)
3707 open (ibond,file=bondname,status='old')
3708 call getenv_loc('THETPAR',thetname)
3709 open (ithep,file=thetname,status='old')
3710 call getenv_loc('ROTPAR',rotname)
3711 open (irotam,file=rotname,status='old')
3712 call getenv_loc('TORPAR',torname)
3713 open (itorp,file=torname,status='old')
3714 call getenv_loc('TORDPAR',tordname)
3715 open (itordp,file=tordname,status='old')
3716 call getenv_loc('SCCORPAR',sccorname)
3717 open (isccor,file=sccorname,status='old')
3718 call getenv_loc('FOURIER',fouriername)
3719 open (ifourier,file=fouriername,status='old')
3720 call getenv_loc('ELEPAR',elename)
3721 open (ielep,file=elename,status='old')
3722 call getenv_loc('SIDEPAR',sidename)
3723 open (isidep,file=sidename,status='old')
3725 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3727 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3728 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3729 ! Get parameter filenames and open the parameter files.
3730 call getenv_loc('BONDPAR',bondname)
3731 open (ibond,file=bondname,status='old',action='read')
3732 call getenv_loc('THETPAR',thetname)
3733 open (ithep,file=thetname,status='old',action='read')
3734 call getenv_loc('ROTPAR',rotname)
3735 open (irotam,file=rotname,status='old',action='read')
3736 call getenv_loc('TORPAR',torname)
3737 open (itorp,file=torname,status='old',action='read')
3738 call getenv_loc('TORDPAR',tordname)
3739 open (itordp,file=tordname,status='old',action='read')
3740 call getenv_loc('SCCORPAR',sccorname)
3741 open (isccor,file=sccorname,status='old',action='read')
3743 call getenv_loc('THETPARPDB',thetname_pdb)
3744 print *,"thetname_pdb ",thetname_pdb
3745 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3746 print *,ithep_pdb," opened"
3748 call getenv_loc('FOURIER',fouriername)
3749 open (ifourier,file=fouriername,status='old',readonly)
3750 call getenv_loc('ELEPAR',elename)
3751 open (ielep,file=elename,status='old',readonly)
3752 call getenv_loc('SIDEPAR',sidename)
3753 open (isidep,file=sidename,status='old',readonly)
3755 call getenv_loc('ROTPARPDB',rotname_pdb)
3756 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3761 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3762 ! Use -DOLDSCP to use hard-coded constants instead.
3764 call getenv_loc('SCPPAR',scpname)
3765 #if defined(WINIFL) || defined(WINPGI)
3766 open (iscpp,file=scpname,status='old',readonly,shared)
3767 #elif (defined CRAY) || (defined AIX)
3768 open (iscpp,file=scpname,status='old',action='read')
3770 open (iscpp,file=scpname,status='old')
3772 open (iscpp,file=scpname,status='old',action='read')
3775 call getenv_loc('PATTERN',patname)
3776 #if defined(WINIFL) || defined(WINPGI)
3777 open (icbase,file=patname,status='old',readonly,shared)
3778 #elif (defined CRAY) || (defined AIX)
3779 open (icbase,file=patname,status='old',action='read')
3781 open (icbase,file=patname,status='old')
3783 open (icbase,file=patname,status='old',action='read')
3786 ! Open output file only for CG processes
3787 ! print *,"Processor",myrank," fg_rank",fg_rank
3788 if (fg_rank.eq.0) then
3790 if (nodes.eq.1) then
3793 npos = dlog10(dfloat(nodes-1))+1
3795 if (npos.lt.3) npos=3
3796 write (liczba,'(i1)') npos
3797 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3799 write (liczba,form) me
3800 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3801 liczba(:ilen(liczba))
3802 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3804 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3806 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3807 liczba(:ilen(liczba))//'.mol2'
3808 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3809 liczba(:ilen(liczba))//'.stat'
3811 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3812 //liczba(:ilen(liczba))//'.stat')
3813 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3816 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3817 liczba(:ilen(liczba))//'.const'
3822 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3823 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3824 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3825 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3826 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3828 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3830 rest2name=prefix(:ilen(prefix))//'.rst'
3832 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3835 #if defined(AIX) || defined(PGI)
3836 if (me.eq.king .or. .not. out1file) &
3837 open(iout,file=outname,status='unknown')
3839 if (fg_rank.gt.0) then
3840 write (liczba,'(i3.3)') myrank/nfgtasks
3841 write (ll,'(bz,i3.3)') fg_rank
3842 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3847 open(igeom,file=intname,status='unknown',position='append')
3848 open(ipdb,file=pdbname,status='unknown')
3849 open(imol2,file=mol2name,status='unknown')
3850 open(istat,file=statname,status='unknown',position='append')
3852 !1out open(iout,file=outname,status='unknown')
3855 if (me.eq.king .or. .not.out1file) &
3856 open(iout,file=outname,status='unknown')
3858 if (fg_rank.gt.0) then
3859 write (liczba,'(i3.3)') myrank/nfgtasks
3860 write (ll,'(bz,i3.3)') fg_rank
3861 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3866 open(igeom,file=intname,status='unknown',access='append')
3867 open(ipdb,file=pdbname,status='unknown')
3868 open(imol2,file=mol2name,status='unknown')
3869 open(istat,file=statname,status='unknown',access='append')
3871 !1out open(iout,file=outname,status='unknown')
3874 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3875 csa_seed=prefix(:lenpre)//'.CSA.seed'
3876 csa_history=prefix(:lenpre)//'.CSA.history'
3877 csa_bank=prefix(:lenpre)//'.CSA.bank'
3878 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3879 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3880 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3881 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3882 csa_int=prefix(:lenpre)//'.int'
3883 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3884 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3885 csa_in=prefix(:lenpre)//'.CSA.in'
3886 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3889 write (iout,'(80(1h-))')
3890 write (iout,'(30x,a)') "FILE ASSIGNMENT"
3891 write (iout,'(80(1h-))')
3892 write (iout,*) "Input file : ",&
3893 pref_orig(:ilen(pref_orig))//'.inp'
3894 write (iout,*) "Output file : ",&
3895 outname(:ilen(outname))
3897 write (iout,*) "Sidechain potential file : ",&
3898 sidename(:ilen(sidename))
3900 write (iout,*) "SCp potential file : ",&
3901 scpname(:ilen(scpname))
3903 write (iout,*) "Electrostatic potential file : ",&
3904 elename(:ilen(elename))
3905 write (iout,*) "Cumulant coefficient file : ",&
3906 fouriername(:ilen(fouriername))
3907 write (iout,*) "Torsional parameter file : ",&
3908 torname(:ilen(torname))
3909 write (iout,*) "Double torsional parameter file : ",&
3910 tordname(:ilen(tordname))
3911 write (iout,*) "SCCOR parameter file : ",&
3912 sccorname(:ilen(sccorname))
3913 write (iout,*) "Bond & inertia constant file : ",&
3914 bondname(:ilen(bondname))
3915 write (iout,*) "Bending parameter file : ",&
3916 thetname(:ilen(thetname))
3917 write (iout,*) "Rotamer parameter file : ",&
3918 rotname(:ilen(rotname))
3921 write (iout,*) "Thetpdb parameter file : ",&
3922 thetname_pdb(:ilen(thetname_pdb))
3925 write (iout,*) "Threading database : ",&
3926 patname(:ilen(patname))
3928 write (iout,*)" DIRTMP : ",&
3930 write (iout,'(80(1h-))')
3933 end subroutine openunits
3934 !-----------------------------------------------------------------------------
3937 use geometry_data, only: nres,dc
3938 use energy_data, only: usampl,iset
3940 ! implicit real*8 (a-h,o-z)
3941 ! include 'DIMENSIONS'
3942 ! include 'COMMON.CHAIN'
3943 ! include 'COMMON.IOUNITS'
3944 ! include 'COMMON.MD'
3947 ! real(kind=8) :: var,ene
3949 open(irest2,file=rest2name,status='unknown')
3950 read(irest2,*) totT,EK,potE,totE,t_bath
3952 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
3955 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
3958 read (irest2,*) iset
3962 end subroutine readrst
3963 !-----------------------------------------------------------------------------
3964 subroutine read_fragments
3968 use control_data, only:out1file
3971 ! implicit real*8 (a-h,o-z)
3972 ! include 'DIMENSIONS'
3976 ! include 'COMMON.SETUP'
3977 ! include 'COMMON.CHAIN'
3978 ! include 'COMMON.IOUNITS'
3979 ! include 'COMMON.MD'
3980 ! include 'COMMON.CONTROL'
3983 ! real(kind=8) :: var,ene
3985 read(inp,*) nset,nfrag,npair,nfrag_back
3987 !el from module energy
3988 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
3989 if(.not.allocated(wfrag_back)) then
3990 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3991 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3993 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
3994 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
3996 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
3997 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4000 if(me.eq.king.or..not.out1file) &
4001 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4002 " nfrag_back",nfrag_back
4004 read(inp,*) mset(iset)
4006 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4008 if(me.eq.king.or..not.out1file) &
4009 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4010 ifrag(2,i,iset), qinfrag(i,iset)
4013 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4015 if(me.eq.king.or..not.out1file) &
4016 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4017 ipair(2,i,iset), qinpair(i,iset)
4020 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4021 wfrag_back(3,i,iset),&
4022 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4023 if(me.eq.king.or..not.out1file) &
4024 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4025 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4029 end subroutine read_fragments
4030 !-----------------------------------------------------------------------------
4032 !-----------------------------------------------------------------------------
4036 ! implicit real*8 (a-h,o-z)
4037 ! include 'DIMENSIONS'
4038 ! include 'COMMON.CSA'
4039 ! include 'COMMON.BANK'
4040 ! include 'COMMON.IOUNITS'
4042 ! integer :: ntf,ik,iw_pdb
4043 ! real(kind=8) :: var,ene
4045 open(icsa_in,file=csa_in,status="old",err=100)
4046 read(icsa_in,*) nconf
4047 read(icsa_in,*) jstart,jend
4048 read(icsa_in,*) nstmax
4049 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4050 read(icsa_in,*) nran0,nran1,irr
4051 read(icsa_in,*) nseed
4052 read(icsa_in,*) ntotal,cut1,cut2
4053 read(icsa_in,*) estop
4054 read(icsa_in,*) icmax,irestart
4055 read(icsa_in,*) ntbankm,dele,difcut
4056 read(icsa_in,*) iref,rmscut,pnccut
4057 read(icsa_in,*) ndiff
4064 end subroutine csa_read
4065 !-----------------------------------------------------------------------------
4066 subroutine initial_write
4069 ! implicit real*8 (a-h,o-z)
4070 ! include 'DIMENSIONS'
4071 ! include 'COMMON.CSA'
4072 ! include 'COMMON.BANK'
4073 ! include 'COMMON.IOUNITS'
4075 ! integer :: ntf,ik,iw_pdb
4076 ! real(kind=8) :: var,ene
4078 open(icsa_seed,file=csa_seed,status="unknown")
4079 write(icsa_seed,*) "seed"
4081 #if defined(AIX) || defined(PGI)
4082 open(icsa_history,file=csa_history,status="unknown",&
4085 open(icsa_history,file=csa_history,status="unknown",&
4088 write(icsa_history,*) nconf
4089 write(icsa_history,*) jstart,jend
4090 write(icsa_history,*) nstmax
4091 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4092 write(icsa_history,*) nran0,nran1,irr
4093 write(icsa_history,*) nseed
4094 write(icsa_history,*) ntotal,cut1,cut2
4095 write(icsa_history,*) estop
4096 write(icsa_history,*) icmax,irestart
4097 write(icsa_history,*) ntbankm,dele,difcut
4098 write(icsa_history,*) iref,rmscut,pnccut
4099 write(icsa_history,*) ndiff
4101 write(icsa_history,*)
4104 open(icsa_bank1,file=csa_bank1,status="unknown")
4105 write(icsa_bank1,*) 0
4109 end subroutine initial_write
4110 !-----------------------------------------------------------------------------
4111 subroutine restart_write
4114 ! implicit real*8 (a-h,o-z)
4115 ! include 'DIMENSIONS'
4116 ! include 'COMMON.IOUNITS'
4117 ! include 'COMMON.CSA'
4118 ! include 'COMMON.BANK'
4120 ! integer :: ntf,ik,iw_pdb
4121 ! real(kind=8) :: var,ene
4123 #if defined(AIX) || defined(PGI)
4124 open(icsa_history,file=csa_history,position="append")
4126 open(icsa_history,file=csa_history,access="append")
4128 write(icsa_history,*)
4129 write(icsa_history,*) "This is restart"
4130 write(icsa_history,*)
4131 write(icsa_history,*) nconf
4132 write(icsa_history,*) jstart,jend
4133 write(icsa_history,*) nstmax
4134 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4135 write(icsa_history,*) nran0,nran1,irr
4136 write(icsa_history,*) nseed
4137 write(icsa_history,*) ntotal,cut1,cut2
4138 write(icsa_history,*) estop
4139 write(icsa_history,*) icmax,irestart
4140 write(icsa_history,*) ntbankm,dele,difcut
4141 write(icsa_history,*) iref,rmscut,pnccut
4142 write(icsa_history,*) ndiff
4143 write(icsa_history,*)
4144 write(icsa_history,*) "irestart is: ", irestart
4146 write(icsa_history,*)
4150 end subroutine restart_write
4151 !-----------------------------------------------------------------------------
4153 !-----------------------------------------------------------------------------
4154 subroutine write_pdb(npdb,titelloc,ee)
4156 ! implicit real*8 (a-h,o-z)
4157 ! include 'DIMENSIONS'
4158 ! include 'COMMON.IOUNITS'
4159 character(len=50) :: titelloc1
4160 character*(*) :: titelloc
4161 character(len=3) :: zahl
4162 character(len=5) :: liczba5
4164 integer :: npdb !,ilen
4168 ! real(kind=8) :: var,ene
4172 if (npdb.lt.1000) then
4173 call numstr(npdb,zahl)
4174 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4176 if (npdb.lt.10000) then
4177 write(liczba5,'(i1,i4)') 0,npdb
4179 write(liczba5,'(i5)') npdb
4181 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4183 call pdbout(ee,titelloc1,ipdb)
4186 end subroutine write_pdb
4187 !-----------------------------------------------------------------------------
4189 !-----------------------------------------------------------------------------
4190 subroutine write_thread_summary
4191 ! Thread the sequence through a database of known structures
4192 use control_data, only: refstr
4194 use energy_data, only: n_ene_comp
4196 ! implicit real*8 (a-h,o-z)
4197 ! include 'DIMENSIONS'
4199 use MPI_data !include 'COMMON.INFO'
4202 ! include 'COMMON.CONTROL'
4203 ! include 'COMMON.CHAIN'
4204 ! include 'COMMON.DBASE'
4205 ! include 'COMMON.INTERACT'
4206 ! include 'COMMON.VAR'
4207 ! include 'COMMON.THREAD'
4208 ! include 'COMMON.FFIELD'
4209 ! include 'COMMON.SBRIDGE'
4210 ! include 'COMMON.HEADER'
4211 ! include 'COMMON.NAMES'
4212 ! include 'COMMON.IOUNITS'
4213 ! include 'COMMON.TIME1'
4215 integer,dimension(maxthread) :: ip
4216 real(kind=8),dimension(0:n_ene) :: energia
4218 integer :: i,j,ii,jj,ipj,ik,kk,ist
4219 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4221 write (iout,'(30x,a/)') &
4222 ' *********** Summary threading statistics ************'
4223 write (iout,'(a)') 'Initial energies:'
4224 write (iout,'(a4,2x,a12,14a14,3a8)') &
4225 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4226 'RMSnat','NatCONT','NNCONT','RMS'
4227 ! Energy sort patterns
4232 enet=ener(n_ene-1,ip(i))
4235 if (ener(n_ene-1,ip(j)).lt.enet) then
4237 enet=ener(n_ene-1,ip(j))
4249 ist=nres_base(2,ii)+ipatt(2,i)
4251 energia(i)=ener0(kk,i)
4253 etot=ener0(n_ene_comp+1,i)
4254 rmsnat=ener0(n_ene_comp+2,i)
4255 rms=ener0(n_ene_comp+3,i)
4256 frac=ener0(n_ene_comp+4,i)
4257 frac_nn=ener0(n_ene_comp+5,i)
4260 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4261 i,str_nam(ii),ist+1,&
4262 (energia(print_order(kk)),kk=1,nprint_ene),&
4263 etot,rmsnat,frac,frac_nn,rms
4265 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4266 i,str_nam(ii),ist+1,&
4267 (energia(print_order(kk)),kk=1,nprint_ene),etot
4270 write (iout,'(//a)') 'Final energies:'
4271 write (iout,'(a4,2x,a12,17a14,3a8)') &
4272 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4273 'RMSnat','NatCONT','NNCONT','RMS'
4277 ist=nres_base(2,ii)+ipatt(2,i)
4279 energia(kk)=ener(kk,ik)
4281 etot=ener(n_ene_comp+1,i)
4282 rmsnat=ener(n_ene_comp+2,i)
4283 rms=ener(n_ene_comp+3,i)
4284 frac=ener(n_ene_comp+4,i)
4285 frac_nn=ener(n_ene_comp+5,i)
4286 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4287 i,str_nam(ii),ist+1,&
4288 (energia(print_order(kk)),kk=1,nprint_ene),&
4289 etot,rmsnat,frac,frac_nn,rms
4291 write (iout,'(/a/)') 'IEXAM array:'
4292 write (iout,'(i5)') nexcl
4294 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4296 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4297 'Max. time for threading step ',max_time_for_thread,&
4298 'Average time for threading step: ',ave_time_for_thread
4300 end subroutine write_thread_summary
4301 !-----------------------------------------------------------------------------
4302 subroutine write_stat_thread(ithread,ipattern,ist)
4304 use energy_data, only: n_ene_comp
4306 ! implicit real*8 (a-h,o-z)
4307 ! include "DIMENSIONS"
4308 ! include "COMMON.CONTROL"
4309 ! include "COMMON.IOUNITS"
4310 ! include "COMMON.THREAD"
4311 ! include "COMMON.FFIELD"
4312 ! include "COMMON.DBASE"
4313 ! include "COMMON.NAMES"
4314 real(kind=8),dimension(0:n_ene) :: energia
4316 integer :: ithread,ipattern,ist,i
4317 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4319 #if defined(AIX) || defined(PGI)
4320 open(istat,file=statname,position='append')
4322 open(istat,file=statname,access='append')
4325 energia(i)=ener(i,ithread)
4327 etot=ener(n_ene_comp+1,ithread)
4328 rmsnat=ener(n_ene_comp+2,ithread)
4329 rms=ener(n_ene_comp+3,ithread)
4330 frac=ener(n_ene_comp+4,ithread)
4331 frac_nn=ener(n_ene_comp+5,ithread)
4332 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4333 ithread,str_nam(ipattern),ist+1,&
4334 (energia(print_order(i)),i=1,nprint_ene),&
4335 etot,rmsnat,frac,frac_nn,rms
4338 end subroutine write_stat_thread
4339 !-----------------------------------------------------------------------------
4341 !-----------------------------------------------------------------------------
4342 end module io_config