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,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm(ntyp)) !(ntyp)
780 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
782 allocate(msc(ntyp+1)) !(ntyp+1)
783 allocate(isc(ntyp+1)) !(ntyp+1)
784 allocate(restok(ntyp+1)) !(ntyp+1)
785 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(long_r_sidechain(ntyp))
787 allocate(short_r_sidechain(ntyp))
792 read (ibond,*) vbldp0,akp,mp,ip,pstok
795 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
796 dsc(i) = vbldsc0(1,i)
800 dsc_inv(i)=1.0D0/dsc(i)
804 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
806 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
807 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
808 dsc(i) = vbldsc0(1,i)
812 dsc_inv(i)=1.0D0/dsc(i)
817 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
818 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
820 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
822 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
823 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
825 write (iout,'(13x,3f10.5)') &
826 vbldsc0(j,i),aksc(j,i),abond0(j,i)
830 !----------------------------------------------------
831 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
832 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
833 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
834 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
835 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
836 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
846 allocate(liptranene(ntyp))
847 !C reading lipid parameters
848 write (iout,*) "iliptranpar",iliptranpar
850 read(iliptranpar,*) pepliptran
853 read(iliptranpar,*) liptranene(i)
854 print *,liptranene(i)
860 ! Read the parameters of the probability distribution/energy expression
861 ! of the virtual-bond valence angles theta
864 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
865 (bthet(j,i,1,1),j=1,2)
866 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
867 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
868 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
872 athet(1,i,1,-1)=athet(1,i,1,1)
873 athet(2,i,1,-1)=athet(2,i,1,1)
874 bthet(1,i,1,-1)=-bthet(1,i,1,1)
875 bthet(2,i,1,-1)=-bthet(2,i,1,1)
876 athet(1,i,-1,1)=-athet(1,i,1,1)
877 athet(2,i,-1,1)=-athet(2,i,1,1)
878 bthet(1,i,-1,1)=bthet(1,i,1,1)
879 bthet(2,i,-1,1)=bthet(2,i,1,1)
883 athet(1,i,-1,-1)=athet(1,-i,1,1)
884 athet(2,i,-1,-1)=-athet(2,-i,1,1)
885 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
886 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
887 athet(1,i,-1,1)=athet(1,-i,1,1)
888 athet(2,i,-1,1)=-athet(2,-i,1,1)
889 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
890 bthet(2,i,-1,1)=bthet(2,-i,1,1)
891 athet(1,i,1,-1)=-athet(1,-i,1,1)
892 athet(2,i,1,-1)=athet(2,-i,1,1)
893 bthet(1,i,1,-1)=bthet(1,-i,1,1)
894 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
899 polthet(j,i)=polthet(j,-i)
902 gthet(j,i)=gthet(j,-i)
910 'Parameters of the virtual-bond valence angles:'
911 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
912 ' ATHETA0 ',' A1 ',' A2 ',&
915 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
916 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
918 write (iout,'(/a/9x,5a/79(1h-))') &
919 'Parameters of the expression for sigma(theta_c):',&
920 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
921 ' ALPH3 ',' SIGMA0C '
923 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
924 (polthet(j,i),j=0,3),sigc0(i)
926 write (iout,'(/a/9x,5a/79(1h-))') &
927 'Parameters of the second gaussian:',&
928 ' THETA0 ',' SIGMA0 ',' G1 ',&
931 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
932 sig0(i),(gthet(j,i),j=1,3)
936 'Parameters of the virtual-bond valence angles:'
937 write (iout,'(/a/9x,5a/79(1h-))') &
938 'Coefficients of expansion',&
939 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
940 ' b1*10^1 ',' b2*10^1 '
942 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
943 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
944 (10*bthet(j,i,1,1),j=1,2)
946 write (iout,'(/a/9x,5a/79(1h-))') &
947 'Parameters of the expression for sigma(theta_c):',&
948 ' alpha0 ',' alph1 ',' alph2 ',&
949 ' alhp3 ',' sigma0c '
951 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
952 (polthet(j,i),j=0,3),sigc0(i)
954 write (iout,'(/a/9x,5a/79(1h-))') &
955 'Parameters of the second gaussian:',&
956 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
959 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
960 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
966 ! Read the parameters of Utheta determined from ab initio surfaces
967 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
969 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
970 ntheterm3,nsingle,ndouble
971 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
973 !----------------------------------------------------
974 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
975 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
976 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
977 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
978 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
979 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
980 !(maxtheterm,-maxthetyp1:maxthetyp1,&
981 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
982 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
983 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
984 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
985 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
986 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
987 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
988 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
989 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
990 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
991 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
992 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
993 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
994 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
995 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
996 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
997 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
999 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1001 ithetyp(i)=-ithetyp(-i)
1004 aa0thet(:,:,:,:)=0.0d0
1005 aathet(:,:,:,:,:)=0.0d0
1006 bbthet(:,:,:,:,:,:)=0.0d0
1007 ccthet(:,:,:,:,:,:)=0.0d0
1008 ddthet(:,:,:,:,:,:)=0.0d0
1009 eethet(:,:,:,:,:,:)=0.0d0
1010 ffthet(:,:,:,:,:,:,:)=0.0d0
1011 ggthet(:,:,:,:,:,:,:)=0.0d0
1013 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1015 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1016 ! VAR:1=non-glicyne non-proline 2=proline
1017 ! VAR:negative values for D-aminoacid
1019 do j=-nthetyp,nthetyp
1020 do k=-nthetyp,nthetyp
1021 read (ithep,'(6a)',end=111,err=111) res1
1022 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1023 ! VAR: aa0thet is variable describing the average value of Foureir
1024 ! VAR: expansion series
1025 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1026 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1027 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1028 read (ithep,*,end=111,err=111) &
1029 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1030 read (ithep,*,end=111,err=111) &
1031 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1032 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1033 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1034 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1036 read (ithep,*,end=111,err=111) &
1037 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1038 ffthet(lll,llll,ll,i,j,k,iblock),&
1039 ggthet(llll,lll,ll,i,j,k,iblock),&
1040 ggthet(lll,llll,ll,i,j,k,iblock),&
1041 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1046 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1047 ! coefficients of theta-and-gamma-dependent terms are zero.
1048 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1049 ! RECOMENTDED AFTER VERSION 3.3)
1053 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1054 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1056 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1057 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1060 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1062 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1065 ! AND COMMENT THE LOOPS BELOW
1069 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1070 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1072 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1073 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1076 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1078 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1083 ! Substitution for D aminoacids from symmetry.
1086 do j=-nthetyp,nthetyp
1087 do k=-nthetyp,nthetyp
1088 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1090 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1094 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1095 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1096 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1097 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1103 ffthet(llll,lll,ll,i,j,k,iblock)= &
1104 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1105 ffthet(lll,llll,ll,i,j,k,iblock)= &
1106 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1107 ggthet(llll,lll,ll,i,j,k,iblock)= &
1108 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1109 ggthet(lll,llll,ll,i,j,k,iblock)= &
1110 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1119 ! Control printout of the coefficients of virtual-bond-angle potentials
1122 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1127 write (iout,'(//4a)') &
1128 'Type ',onelett(i),onelett(j),onelett(k)
1129 write (iout,'(//a,10x,a)') " l","a[l]"
1130 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1131 write (iout,'(i2,1pe15.5)') &
1132 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1134 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1135 "b",l,"c",l,"d",l,"e",l
1137 write (iout,'(i2,4(1pe15.5))') m,&
1138 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1139 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1143 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1144 "f+",l,"f-",l,"g+",l,"g-",l
1147 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1148 ffthet(n,m,l,i,j,k,iblock),&
1149 ffthet(m,n,l,i,j,k,iblock),&
1150 ggthet(n,m,l,i,j,k,iblock),&
1151 ggthet(m,n,l,i,j,k,iblock)
1161 write (2,*) "Start reading THETA_PDB",ithep_pdb
1163 ! write (2,*) 'i=',i
1164 read (ithep_pdb,*,err=111,end=111) &
1165 a0thet(i),(athet(j,i,1,1),j=1,2),&
1166 (bthet(j,i,1,1),j=1,2)
1167 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1168 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1169 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1170 sigc0(i)=sigc0(i)**2
1173 athet(1,i,1,-1)=athet(1,i,1,1)
1174 athet(2,i,1,-1)=athet(2,i,1,1)
1175 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1176 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1177 athet(1,i,-1,1)=-athet(1,i,1,1)
1178 athet(2,i,-1,1)=-athet(2,i,1,1)
1179 bthet(1,i,-1,1)=bthet(1,i,1,1)
1180 bthet(2,i,-1,1)=bthet(2,i,1,1)
1183 a0thet(i)=a0thet(-i)
1184 athet(1,i,-1,-1)=athet(1,-i,1,1)
1185 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1186 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1187 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1188 athet(1,i,-1,1)=athet(1,-i,1,1)
1189 athet(2,i,-1,1)=-athet(2,-i,1,1)
1190 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1191 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1192 athet(1,i,1,-1)=-athet(1,-i,1,1)
1193 athet(2,i,1,-1)=athet(2,-i,1,1)
1194 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1195 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1196 theta0(i)=theta0(-i)
1200 polthet(j,i)=polthet(j,-i)
1203 gthet(j,i)=gthet(j,-i)
1206 write (2,*) "End reading THETA_PDB"
1211 !-------------------------------------------
1212 allocate(nlob(ntyp1)) !(ntyp1)
1213 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1214 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1215 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1222 gaussc(:,:,:,:)=0.0D0
1226 ! Read the parameters of the probability distribution/energy expression
1227 ! of the side chains.
1230 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1234 dsc_inv(i)=1.0D0/dsc(i)
1245 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1246 ((blower(k,l,1),l=1,k),k=1,3)
1247 censc(1,1,-i)=censc(1,1,i)
1248 censc(2,1,-i)=censc(2,1,i)
1249 censc(3,1,-i)=-censc(3,1,i)
1251 read (irotam,*,end=112,err=112) bsc(j,i)
1252 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1253 ((blower(k,l,j),l=1,k),k=1,3)
1254 censc(1,j,-i)=censc(1,j,i)
1255 censc(2,j,-i)=censc(2,j,i)
1256 censc(3,j,-i)=-censc(3,j,i)
1257 ! BSC is amplitude of Gaussian
1264 akl=akl+blower(k,m,j)*blower(l,m,j)
1268 if (((k.eq.3).and.(l.ne.3)) &
1269 .or.((l.eq.3).and.(k.ne.3))) then
1270 gaussc(k,l,j,-i)=-akl
1271 gaussc(l,k,j,-i)=-akl
1273 gaussc(k,l,j,-i)=akl
1274 gaussc(l,k,j,-i)=akl
1283 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1286 if (nlobi.gt.0) then
1288 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1289 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1290 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1291 'log h',(bsc(j,i),j=1,nlobi)
1292 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1293 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1295 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1296 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1299 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1300 write (iout,'(a,f10.4,4(16x,f10.4))') &
1301 'Center ',(bsc(j,i),j=1,nlobi)
1302 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1311 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1312 ! added by Urszula Kozlowska 07/11/2007
1314 !el Maximum number of SC local term fitting function coefficiants
1315 !el integer,parameter :: maxsccoef=65
1317 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1320 read (irotam,*,end=112,err=112)
1322 read (irotam,*,end=112,err=112)
1325 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1330 ! Read the parameters of the probability distribution/energy expression
1331 ! of the side chains.
1333 write (2,*) "Start reading ROTAM_PDB"
1335 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1339 dsc_inv(i)=1.0D0/dsc(i)
1350 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1351 ((blower(k,l,1),l=1,k),k=1,3)
1353 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1354 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1355 ((blower(k,l,j),l=1,k),k=1,3)
1362 akl=akl+blower(k,m,j)*blower(l,m,j)
1372 write (2,*) "End reading ROTAM_PDB"
1378 ! Read torsional parameters in old format
1380 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1382 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1383 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1384 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1386 !el from energy module--------
1387 allocate(v1(nterm_old,ntortyp,ntortyp))
1388 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1389 !el---------------------------
1394 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1400 write (iout,'(/a/)') 'Torsional constants:'
1403 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1404 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1410 ! Read torsional parameters
1412 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1413 read (itorp,*,end=113,err=113) ntortyp
1414 !el from energy module---------
1415 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1416 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1418 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1419 allocate(vlor2(maxlor,ntortyp,ntortyp))
1420 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1421 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1423 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1424 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1425 !el---------------------------
1428 !el---------------------------
1430 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1432 itortyp(i)=-itortyp(-i)
1434 itortyp(ntyp1)=ntortyp
1435 itortyp(-ntyp1)=-ntortyp
1437 write (iout,*) 'ntortyp',ntortyp
1439 do j=-ntortyp+1,ntortyp-1
1440 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1442 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1443 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1446 do k=1,nterm(i,j,iblock)
1447 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1449 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1450 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1451 v0ij=v0ij+si*v1(k,i,j,iblock)
1453 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1454 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1455 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1457 do k=1,nlor(i,j,iblock)
1458 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1459 vlor2(k,i,j),vlor3(k,i,j)
1460 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1463 v0(-i,-j,iblock)=v0ij
1469 write (iout,'(/a/)') 'Torsional constants:'
1471 do i=-ntortyp,ntortyp
1472 do j=-ntortyp,ntortyp
1473 write (iout,*) 'ityp',i,' jtyp',j
1474 write (iout,*) 'Fourier constants'
1475 do k=1,nterm(i,j,iblock)
1476 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1479 write (iout,*) 'Lorenz constants'
1480 do k=1,nlor(i,j,iblock)
1481 write (iout,'(3(1pe15.5))') &
1482 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1488 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1490 ! 6/23/01 Read parameters for double torsionals
1492 !el from energy module------------
1493 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1494 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1495 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1496 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1497 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1498 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1499 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1500 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1501 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1502 !---------------------------------
1506 do j=-ntortyp+1,ntortyp-1
1507 do k=-ntortyp+1,ntortyp-1
1508 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1509 ! write (iout,*) "OK onelett",
1512 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1513 .or. t3.ne.toronelet(k)) then
1514 write (iout,*) "Error in double torsional parameter file",&
1517 call MPI_Finalize(Ierror)
1519 stop "Error in double torsional parameter file"
1521 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1522 ntermd_2(i,j,k,iblock)
1523 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1524 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1525 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1526 ntermd_1(i,j,k,iblock))
1527 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1528 ntermd_1(i,j,k,iblock))
1529 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1530 ntermd_1(i,j,k,iblock))
1531 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1532 ntermd_1(i,j,k,iblock))
1533 ! Martix of D parameters for one dimesional foureir series
1534 do l=1,ntermd_1(i,j,k,iblock)
1535 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1536 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1537 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1538 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1539 ! write(iout,*) "whcodze" ,
1540 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1542 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1543 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1544 v2s(m,l,i,j,k,iblock),&
1545 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1546 ! Martix of D parameters for two dimesional fourier series
1547 do l=1,ntermd_2(i,j,k,iblock)
1549 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1550 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1551 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1552 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1561 write (iout,*) 'Constants for double torsionals'
1564 do j=-ntortyp+1,ntortyp-1
1565 do k=-ntortyp+1,ntortyp-1
1566 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1567 ' nsingle',ntermd_1(i,j,k,iblock),&
1568 ' ndouble',ntermd_2(i,j,k,iblock)
1570 write (iout,*) 'Single angles:'
1571 do l=1,ntermd_1(i,j,k,iblock)
1572 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1573 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1574 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1575 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1578 write (iout,*) 'Pairs of angles:'
1579 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1580 do l=1,ntermd_2(i,j,k,iblock)
1581 write (iout,'(i5,20f10.5)') &
1582 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1585 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1586 do l=1,ntermd_2(i,j,k,iblock)
1587 write (iout,'(i5,20f10.5)') &
1588 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1589 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1598 ! Read of Side-chain backbone correlation parameters
1599 ! Modified 11 May 2012 by Adasko
1602 read (isccor,*,end=119,err=119) nsccortyp
1604 !el from module energy-------------
1605 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1606 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1607 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1608 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1609 !-----------------------------------
1611 !el from module energy-------------
1612 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1614 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1616 isccortyp(i)=-isccortyp(-i)
1618 iscprol=isccortyp(20)
1619 ! write (iout,*) 'ntortyp',ntortyp
1621 !c maxinter is maximum interaction sites
1622 !el from module energy---------
1623 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1624 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1625 -nsccortyp:nsccortyp))
1626 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1627 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1628 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1629 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1630 !-----------------------------------
1634 read (isccor,*,end=119,err=119) &
1635 nterm_sccor(i,j),nlor_sccor(i,j)
1641 nterm_sccor(-i,j)=nterm_sccor(i,j)
1642 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1643 nterm_sccor(i,-j)=nterm_sccor(i,j)
1644 do k=1,nterm_sccor(i,j)
1645 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1647 if (j.eq.iscprol) then
1648 if (i.eq.isccortyp(10)) then
1649 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1650 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1652 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1653 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1654 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1655 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1656 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1657 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1658 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1659 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1662 if (i.eq.isccortyp(10)) then
1663 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1664 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1666 if (j.eq.isccortyp(10)) then
1667 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1668 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1670 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1671 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1672 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1673 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1674 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1675 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1679 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1680 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1681 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1682 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1685 do k=1,nlor_sccor(i,j)
1686 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1687 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1688 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1689 (1+vlor3sccor(k,i,j)**2)
1691 v0sccor(l,i,j)=v0ijsccor
1692 v0sccor(l,-i,j)=v0ijsccor1
1693 v0sccor(l,i,-j)=v0ijsccor2
1694 v0sccor(l,-i,-j)=v0ijsccor3
1700 !el from module energy-------------
1701 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1703 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1704 ! write (iout,*) 'ntortyp',ntortyp
1706 !c maxinter is maximum interaction sites
1707 !el from module energy---------
1708 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1709 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1710 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1711 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1712 !-----------------------------------
1716 read (isccor,*,end=119,err=119) &
1717 nterm_sccor(i,j),nlor_sccor(i,j)
1721 do k=1,nterm_sccor(i,j)
1722 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1724 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1727 do k=1,nlor_sccor(i,j)
1728 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1729 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1730 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1731 (1+vlor3sccor(k,i,j)**2)
1733 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1741 write (iout,'(/a/)') 'Torsional constants:'
1744 write (iout,*) 'ityp',i,' jtyp',j
1745 write (iout,*) 'Fourier constants'
1746 do k=1,nterm_sccor(i,j)
1747 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1749 write (iout,*) 'Lorenz constants'
1750 do k=1,nlor_sccor(i,j)
1751 write (iout,'(3(1pe15.5))') &
1752 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1759 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1760 ! interaction energy of the Gly, Ala, and Pro prototypes.
1765 write (iout,*) "Coefficients of the cumulants"
1767 read (ifourier,*) nloctyp
1768 !write(iout,*) "nloctyp",nloctyp
1769 !el from module energy-------
1770 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1771 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1772 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1773 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1774 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1775 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1776 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1777 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1781 !--------------------------------
1784 read (ifourier,*,end=115,err=115)
1785 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1787 write (iout,*) 'Type',i
1788 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1798 B1tilde(1,-i) =-b(3)
1800 ! b1tilde(1,i)=0.0d0
1801 ! b1tilde(2,i)=0.0d0
1826 Ctilde(1,2,-i)=-b(9)
1830 ! Ctilde(1,1,i)=0.0d0
1831 ! Ctilde(1,2,i)=0.0d0
1832 ! Ctilde(2,1,i)=0.0d0
1833 ! Ctilde(2,2,i)=0.0d0
1851 Dtilde(1,2,-i)=-b(8)
1855 ! Dtilde(1,1,i)=0.0d0
1856 ! Dtilde(1,2,i)=0.0d0
1857 ! Dtilde(2,1,i)=0.0d0
1858 ! Dtilde(2,2,i)=0.0d0
1859 EE(1,1,i)= b(10)+b(11)
1860 EE(2,2,i)=-b(10)+b(11)
1861 EE(2,1,i)= b(12)-b(13)
1862 EE(1,2,i)= b(12)+b(13)
1863 EE(1,1,-i)= b(10)+b(11)
1864 EE(2,2,-i)=-b(10)+b(11)
1865 EE(2,1,-i)=-b(12)+b(13)
1866 EE(1,2,-i)=-b(12)-b(13)
1872 ! ee(2,1,i)=ee(1,2,i)
1876 write (iout,*) 'Type',i
1878 write(iout,*) B1(1,i),B1(2,i)
1880 write(iout,*) B2(1,i),B2(2,i)
1883 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1887 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1891 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1896 ! Read electrostatic-interaction parameters
1901 write (iout,'(/a)') 'Electrostatic interaction constants:'
1902 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1903 'IT','JT','APP','BPP','AEL6','AEL3'
1905 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1906 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1907 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1908 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1913 app (i,j)=epp(i,j)*rri*rri
1914 bpp (i,j)=-2.0D0*epp(i,j)*rri
1915 ael6(i,j)=elpp6(i,j)*4.2D0**6
1916 ael3(i,j)=elpp3(i,j)*4.2D0**3
1918 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1924 ! Read side-chain interaction parameters.
1926 !el from module energy - COMMON.INTERACT-------
1927 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1928 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1929 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1930 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1931 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1932 allocate(epslip(ntyp,ntyp))
1940 !--------------------------------
1942 read (isidep,*,end=117,err=117) ipot,expon
1943 if (ipot.lt.1 .or. ipot.gt.5) then
1944 write (iout,'(2a)') 'Error while reading SC interaction',&
1945 'potential file - unknown potential type.'
1947 call MPI_Finalize(Ierror)
1953 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1954 ', exponents are ',expon,2*expon
1955 ! goto (10,20,30,30,40) ipot
1957 !----------------------- LJ potential ---------------------------------
1959 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1960 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961 (sigma0(i),i=1,ntyp)
1963 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1964 write (iout,'(a/)') 'The epsilon array:'
1965 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1966 write (iout,'(/a)') 'One-body parameters:'
1967 write (iout,'(a,4x,a)') 'residue','sigma'
1968 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
1971 !----------------------- LJK potential --------------------------------
1973 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1974 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1975 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1977 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1978 write (iout,'(a/)') 'The epsilon array:'
1979 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1980 write (iout,'(/a)') 'One-body parameters:'
1981 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1982 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
1986 !---------------------- GB or BP potential -----------------------------
1990 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1992 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1993 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1994 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1995 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1997 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2000 ! For the GB potential convert sigma'**2 into chi'
2003 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2007 write (iout,'(/a/)') 'Parameters of the BP potential:'
2008 write (iout,'(a/)') 'The epsilon array:'
2009 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010 write (iout,'(/a)') 'One-body parameters:'
2011 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2013 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2014 chip(i),alp(i),i=1,ntyp)
2017 !--------------------- GBV potential -----------------------------------
2019 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2020 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2021 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2022 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2024 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2025 write (iout,'(a/)') 'The epsilon array:'
2026 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2027 write (iout,'(/a)') 'One-body parameters:'
2028 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2029 's||/s_|_^2',' chip ',' alph '
2030 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2031 sigii(i),chip(i),alp(i),i=1,ntyp)
2034 write(iout,*)"Wrong ipot"
2039 !-----------------------------------------------------------------------
2040 ! Calculate the "working" parameters of SC interactions.
2042 !el from module energy - COMMON.INTERACT-------
2043 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2044 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2045 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2046 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2048 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2061 sc_aa_tube_par(:)=0.0d0
2062 sc_bb_tube_par(:)=0.0d0
2064 !--------------------------------
2069 epslip(i,j)=epslip(j,i)
2074 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2075 sigma(j,i)=sigma(i,j)
2076 rs0(i,j)=dwa16*sigma(i,j)
2080 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2081 'Working parameters of the SC interactions:',&
2082 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2087 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2096 sigeps=dsign(1.0D0,epsij)
2098 aa_aq(i,j)=epsij*rrij*rrij
2099 bb_aq(i,j)=-sigeps*epsij*rrij
2100 aa_aq(j,i)=aa_aq(i,j)
2101 bb_aq(j,i)=bb_aq(i,j)
2102 epsijlip=epslip(i,j)
2103 sigeps=dsign(1.0D0,epsijlip)
2104 epsijlip=dabs(epsijlip)
2105 aa_lip(i,j)=epsijlip*rrij*rrij
2106 bb_lip(i,j)=-sigeps*epsijlip*rrij
2107 aa_lip(j,i)=aa_lip(i,j)
2108 bb_lip(j,i)=bb_lip(i,j)
2109 !C write(iout,*) aa_lip
2111 sigt1sq=sigma0(i)**2
2112 sigt2sq=sigma0(j)**2
2115 ratsig1=sigt2sq/sigt1sq
2116 ratsig2=1.0D0/ratsig1
2117 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2118 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2119 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2123 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2124 sigmaii(i,j)=rsum_max
2125 sigmaii(j,i)=rsum_max
2127 ! sigmaii(i,j)=r0(i,j)
2128 ! sigmaii(j,i)=r0(i,j)
2130 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2131 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2132 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2133 augm(i,j)=epsij*r_augm**(2*expon)
2134 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2141 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2142 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2143 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2147 write(iout,*) "tube param"
2148 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2149 ccavtubpep,dcavtubpep,tubetranenepep
2150 sigmapeptube=sigmapeptube**6
2151 sigeps=dsign(1.0D0,epspeptube)
2152 epspeptube=dabs(epspeptube)
2153 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2154 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2155 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2157 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2158 ccavtub(i),dcavtub(i),tubetranene(i)
2159 sigmasctube=sigmasctube**6
2160 sigeps=dsign(1.0D0,epssctube)
2161 epssctube=dabs(epssctube)
2162 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2163 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2164 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2167 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2172 ! Define the SC-p interaction constants (hard-coded; old style)
2175 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2177 ! aad(i,1)=0.3D0*4.0D0**12
2178 ! Following line for constants currently implemented
2179 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2180 aad(i,1)=1.5D0*4.0D0**12
2181 ! aad(i,1)=0.17D0*5.6D0**12
2183 ! "Soft" SC-p repulsion
2185 ! Following line for constants currently implemented
2186 ! aad(i,1)=0.3D0*4.0D0**6
2187 ! "Hard" SC-p repulsion
2188 bad(i,1)=3.0D0*4.0D0**6
2189 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2198 ! 8/9/01 Read the SC-p interaction constants from file
2201 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2204 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2205 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2206 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2207 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2211 write (iout,*) "Parameters of SC-p interactions:"
2213 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2214 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2220 ! Define the constants of the disulfide bridge
2224 ! Old arbitrary potential - commented out.
2229 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2230 ! energy surface of diethyl disulfide.
2231 ! A. Liwo and U. Kozlowska, 11/24/03
2248 write (iout,'(/a)') "Disulfide bridge parameters:"
2249 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2250 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2251 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2252 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2256 111 write (iout,*) "Error reading bending energy parameters."
2258 112 write (iout,*) "Error reading rotamer energy parameters."
2260 113 write (iout,*) "Error reading torsional energy parameters."
2262 114 write (iout,*) "Error reading double torsional energy parameters."
2264 115 write (iout,*) &
2265 "Error reading cumulant (multibody energy) parameters."
2267 116 write (iout,*) "Error reading electrostatic energy parameters."
2269 117 write (iout,*) "Error reading side chain interaction parameters."
2271 118 write (iout,*) "Error reading SCp interaction parameters."
2273 119 write (iout,*) "Error reading SCCOR parameters"
2276 call MPI_Finalize(Ierror)
2280 end subroutine parmread
2282 !-----------------------------------------------------------------------------
2284 !-----------------------------------------------------------------------------
2285 subroutine printmat(ldim,m,n,iout,key,a)
2288 character(len=3),dimension(n) :: key
2289 real(kind=8),dimension(ldim,n) :: a
2291 integer :: i,j,k,m,iout,nlim
2295 write (iout,1000) (key(k),k=i,nlim)
2297 1000 format (/5x,8(6x,a3))
2298 1020 format (/80(1h-)/)
2300 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2303 1010 format (a3,2x,8(f9.4))
2305 end subroutine printmat
2306 !-----------------------------------------------------------------------------
2308 !-----------------------------------------------------------------------------
2310 ! Read the PDB file and convert the peptide geometry into virtual-chain
2313 use energy_data, only: itype
2317 use control, only: rescode,sugarcode
2318 ! implicit real*8 (a-h,o-z)
2319 ! include 'DIMENSIONS'
2320 ! include 'COMMON.LOCAL'
2321 ! include 'COMMON.VAR'
2322 ! include 'COMMON.CHAIN'
2323 ! include 'COMMON.INTERACT'
2324 ! include 'COMMON.IOUNITS'
2325 ! include 'COMMON.GEO'
2326 ! include 'COMMON.NAMES'
2327 ! include 'COMMON.CONTROL'
2328 ! include 'COMMON.DISTFIT'
2329 ! include 'COMMON.SETUP'
2330 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2332 logical :: lprn=.true.,fail
2333 real(kind=8),dimension(3) :: e1,e2,e3
2334 real(kind=8) :: dcj,efree_temp
2335 character(len=3) :: seq,res
2336 character(len=5) :: atom
2337 character(len=80) :: card
2338 real(kind=8),dimension(3,20) :: sccor
2339 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2341 character*1 :: sugar
2343 real(kind=8),dimension(3) :: ccc
2345 integer,dimension(2,maxres/3) :: hfrag_alloc
2346 integer,dimension(4,maxres/3) :: bfrag_alloc
2347 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2348 real(kind=8),dimension(:,:), allocatable :: c_temporary
2349 integer,dimension(:,:) , allocatable :: itype_temporary
2356 ! write (2,*) "UNRES_PDB",unres_pdb
2364 !-----------------------------
2365 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2366 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2369 read (ipdbin,'(a80)',end=10) card
2370 ! write (iout,'(a)') card
2371 if (card(:5).eq.'HELIX') then
2374 read(card(22:25),*) hfrag(1,nhfrag)
2375 read(card(34:37),*) hfrag(2,nhfrag)
2377 if (card(:5).eq.'SHEET') then
2380 read(card(24:26),*) bfrag(1,nbfrag)
2381 read(card(35:37),*) bfrag(2,nbfrag)
2382 !rc----------------------------------------
2383 !rc to be corrected !!!
2384 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2385 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2386 !rc----------------------------------------
2388 if (card(:3).eq.'END') then
2390 else if (card(:3).eq.'TER') then
2395 itype(ires_old,molecule)=ntyp1_molec(molecule)
2396 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2397 nres_molec(molecule)=nres_molec(molecule)+2
2399 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2402 dc(j,ires)=sccor(j,iii)
2405 call sccenter(ires,iii,sccor)
2411 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2412 ! Fish out the ATOM cards.
2413 ! write(iout,*) 'card',card(1:20)
2414 if (index(card(1:4),'ATOM').gt.0) then
2415 read (card(12:16),*) atom
2416 ! write (iout,*) "! ",atom," !",ires
2417 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2418 read (card(23:26),*) ires
2419 read (card(18:20),'(a3)') res
2420 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2421 ! & " ires_old",ires_old
2422 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2423 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2424 if (ires-ishift+ishift1.ne.ires_old) then
2425 ! Calculate the CM of the preceding residue.
2426 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2428 ! write (iout,*) "Calculating sidechain center iii",iii
2431 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2434 call sccenter(ires_old,iii,sccor)
2438 ! Start new residue.
2439 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2442 else if (ibeg.eq.1) then
2443 write (iout,*) "BEG ires",ires
2445 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2448 nres_molec(molecule)=nres_molec(molecule)+1
2450 ires=ires-ishift+ishift1
2452 ! write (iout,*) "ishift",ishift," ires",ires,&
2453 ! " ires_old",ires_old
2455 else if (ibeg.eq.2) then
2457 ishift=-ires_old+ires-1 !!!!!
2458 ishift1=ishift1-1 !!!!!
2459 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2460 ires=ires-ishift+ishift1
2461 print *,ires,ishift,ishift1
2465 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2466 ires=ires-ishift+ishift1
2469 ! print *,'atom',ires,atom
2470 if (res.eq.'ACE' .or. res.eq.'NHE') then
2473 if (atom.eq.'CA '.or.atom.eq.'N ') then
2475 itype(ires,molecule)=rescode(ires,res,0,molecule)
2476 ! nres_molec(molecule)=nres_molec(molecule)+1
2479 itype(ires,molecule)=rescode(ires,res(3:4),0,molecule)
2480 ! nres_molec(molecule)=nres_molec(molecule)+1
2484 ires=ires-ishift+ishift1
2486 ! write (iout,*) "ires_old",ires_old," ires",ires
2487 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2490 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2491 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2492 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2493 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2494 ! print *,ires,ishift,ishift1
2495 ! write (iout,*) "backbone ",atom
2497 write (iout,'(2i3,2x,a,3f8.3)') &
2498 ires,itype(ires,1),res,(c(j,ires),j=1,3)
2501 nres_molec(molecule)=nres_molec(molecule)+1
2503 sccor(j,iii)=c(j,ires)
2505 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2506 atom.eq."C2'" .or. atom.eq."C3'" &
2507 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2508 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2509 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2510 print *,ires,ishift,ishift1
2514 ! sccor(j,iii)=c(j,ires)
2517 c(j,ires)=c(j,ires)+ccc(j)/5.0
2519 if (counter.eq.5) then
2521 nres_molec(molecule)=nres_molec(molecule)+1
2523 ! sccor(j,iii)=c(j,ires)
2527 ! print *, "ATOM",atom(1:3)
2528 else if (atom(1:3).eq."C5'") then
2529 read (card(19:19),'(a1)') sugar
2530 isugar=sugarcode(sugar,ires)
2537 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2540 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2542 ! write (*,*) card(23:27),ires,itype(ires,1)
2543 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2544 atom.ne.'N' .and. atom.ne.'C' .and. &
2545 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2546 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2547 .and. atom.ne.'P '.and. &
2548 atom(1:1).ne.'H' .and. &
2549 atom.ne.'OP1' .and. atom.ne.'OP2 ') then
2550 ! write (iout,*) "sidechain ",atom
2551 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2552 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2553 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2555 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2560 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2561 if (ires.eq.0) return
2562 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2565 if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
2566 ! print *,'I have', nres_molec(:)
2569 if (nres_molec(k).eq.0) cycle
2571 ! write (iout,*) i,itype(i,1)
2572 ! if (itype(i,1).eq.ntyp1) then
2573 ! write (iout,*) "dummy",i,itype(i,1)
2575 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2576 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2580 if (itype(i,k).eq.ntyp1_molec(k)) then
2581 if (itype(i+1,k).eq.ntyp1_molec(k)) then
2582 if (itype(i+2,k).eq.0) then
2583 print *,"masz sieczke"
2585 if (itype(i+2,j).ne.0) then
2587 itype(i+1,j)=ntyp1_molec(j)
2588 nres_molec(k)=nres_molec(k)-1
2589 nres_molec(j)=nres_molec(j)+1
2595 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2596 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2597 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2599 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2600 ! print *,i,'tu dochodze'
2601 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2609 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2613 dcj=(c(j,i-2)-c(j,i-3))/2.0
2614 if (dcj.eq.0) dcj=1.23591524223
2619 else !itype(i+1,1).eq.ntyp1
2621 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2622 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2629 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2633 dcj=(c(j,i+3)-c(j,i+2))/2.0
2634 if (dcj.eq.0) dcj=1.23591524223
2639 endif !itype(i+1,1).eq.ntyp1
2640 endif !itype.eq.ntyp1
2644 ! Calculate the CM of the last side chain.
2648 dc(j,ires)=sccor(j,iii)
2651 call sccenter(ires,iii,sccor)
2657 ! print *,"molecule",molecule
2658 if (itype(nres,1).ne.10) then
2660 itype(nres,molecule)=ntyp1_molec(molecule)
2661 nres_molec(molecule)=nres_molec(molecule)+1
2663 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2664 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2671 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2675 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2676 c(j,nres)=c(j,nres-1)+dcj
2677 c(j,2*nres)=c(j,nres)
2681 ! print *,'I have',nres, nres_molec(:)
2683 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2685 if (nres.ne.nres0) then
2686 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2688 stop "Error nres value in WHAM input"
2691 !---------------------------------
2692 !el reallocate tables
2695 ! hfrag_alloc(j,i)=hfrag(j,i)
2698 ! bfrag_alloc(j,i)=bfrag(j,i)
2704 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2705 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2706 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2707 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2711 ! hfrag(j,i)=hfrag_alloc(j,i)
2716 ! bfrag(j,i)=bfrag_alloc(j,i)
2719 !el end reallocate tables
2720 !---------------------------------
2728 c(j,2*nres)=c(j,nres)
2731 if (itype(1,1).eq.ntyp1) then
2735 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2736 call refsys(2,3,4,e1,e2,e3,fail)
2743 c(j,1)=c(j,2)-1.9d0*e2(j)
2747 dcj=(c(j,4)-c(j,3))/2.0
2753 ! First lets assign correct dummy to correct type of chain
2755 if (itype(1,1).eq.ntyp1) then
2756 if (itype(2,1).eq.0) then
2758 if (itype(2,j).ne.0) then
2760 itype(1,j)=ntyp1_molec(j)
2761 nres_molec(1)=nres_molec(1)-1
2762 nres_molec(j)=nres_molec(j)+1
2769 print *,'I have',nres, nres_molec(:)
2771 ! Copy the coordinates to reference coordinates
2777 ! Calculate internal coordinates.
2779 write (iout,'(/a)') &
2780 "Cartesian coordinates of the reference structure"
2781 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2782 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2784 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
2785 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
2786 (c(j,ires+nres),j=1,3)
2789 ! znamy już nres wiec mozna alokowac tablice
2790 ! Calculate internal coordinates.
2791 if(me.eq.king.or..not.out1file)then
2792 write (iout,'(a)') &
2793 "Backbone and SC coordinates as read from the PDB"
2795 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2796 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
2797 (c(j,nres+ires),j=1,3)
2800 ! NOW LETS ROCK! SORTING
2801 allocate(c_temporary(3,2*nres))
2802 allocate(itype_temporary(nres,5))
2803 itype_temporary(:,:)=0
2807 if (itype(i,k).ne.0) then
2809 c_temporary(j,seqalingbegin)=c(j,i)
2810 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
2813 itype_temporary(seqalingbegin,k)=itype(i,k)
2814 seqalingbegin=seqalingbegin+1
2820 c(j,i)=c_temporary(j,i)
2825 itype(i,k)=itype_temporary(i,k)
2828 if (itype(1,1).eq.ntyp1) then
2832 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2833 call refsys(2,3,4,e1,e2,e3,fail)
2840 c(j,1)=c(j,2)-1.9d0*e2(j)
2844 dcj=(c(j,4)-c(j,3))/2.0
2852 write (iout,'(/a)') &
2853 "Cartesian coordinates of the reference structure after sorting"
2854 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2855 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2857 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
2858 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
2859 (c(j,ires+nres),j=1,3)
2863 print *,seqalingbegin,nres
2864 if(.not.allocated(vbld)) then
2865 allocate(vbld(2*nres))
2870 if(.not.allocated(vbld_inv)) then
2871 allocate(vbld_inv(2*nres))
2877 if(.not.allocated(theta)) then
2878 allocate(theta(nres+2))
2882 if(.not.allocated(phi)) allocate(phi(nres+2))
2883 if(.not.allocated(alph)) allocate(alph(nres+2))
2884 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2885 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2886 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2887 if(.not.allocated(costtab)) allocate(costtab(nres))
2888 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2889 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2890 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2891 if(.not.allocated(xxref)) allocate(xxref(nres))
2892 if(.not.allocated(yyref)) allocate(yyref(nres))
2893 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2894 if(.not.allocated(dc_norm)) then
2895 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2896 allocate(dc_norm(3,0:2*nres+2))
2900 call int_from_cart(.true.,.false.)
2901 call sc_loc_geom(.false.)
2903 thetaref(i)=theta(i)
2913 dc(j,i)=c(j,i+1)-c(j,i)
2914 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2919 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2920 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2922 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2926 ! Copy the coordinates to reference coordinates
2927 ! Splits to single chain if occurs
2931 ! cref(j,i,cou)=c(j,i)
2935 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2936 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2937 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2938 !-----------------------------
2944 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2946 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
2949 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2954 cref(j,i,cou)=c(j,i)
2955 cref(j,i+nres,cou)=c(j,i+nres)
2957 chain_rep(j,lll,kkk)=c(j,i)
2958 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2962 write (iout,*) chain_length
2963 if (chain_length.eq.0) chain_length=nres
2965 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2966 chain_rep(j,chain_length+nres,symetr) &
2967 =chain_rep(j,chain_length+nres,1)
2970 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2972 ! do kkk=1,chain_length
2973 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2977 ! makes copy of chains
2978 write (iout,*) "symetr", symetr
2983 if (symetr.gt.1) then
2990 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2996 ! write (iout,*) i,icha
2997 do lll=1,chain_length
2999 if (cou.le.nres) then
3001 kupa=mod(lll,chain_length)
3002 iprzes=(kkk-1)*chain_length+lll
3003 if (kupa.eq.0) kupa=chain_length
3004 ! write (iout,*) "kupa", kupa
3005 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3006 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3013 !-koniec robienia kopii
3016 write (iout,*) "nowa struktura", nperm
3018 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3020 cref(3,i,kkk),cref(1,nres+i,kkk),&
3021 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3023 100 format (//' alpha-carbon coordinates ',&
3024 ' centroid coordinates'/ &
3025 ' ', 6X,'X',11X,'Y',11X,'Z', &
3026 10X,'X',11X,'Y',11X,'Z')
3027 110 format (a,'(',i3,')',6f12.5)
3033 bfrag(i,j)=bfrag(i,j)-ishift
3039 hfrag(i,j)=hfrag(i,j)-ishift
3045 end subroutine readpdb
3046 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3047 !-----------------------------------------------------------------------------
3049 !-----------------------------------------------------------------------------
3050 subroutine read_control
3064 use random, only: random_init
3065 ! implicit real*8 (a-h,o-z)
3066 ! include 'DIMENSIONS'
3068 use prng, only:prng_restart
3070 logical :: OKRandom!, prng_restart
3073 ! include 'COMMON.IOUNITS'
3074 ! include 'COMMON.TIME1'
3075 ! include 'COMMON.THREAD'
3076 ! include 'COMMON.SBRIDGE'
3077 ! include 'COMMON.CONTROL'
3078 ! include 'COMMON.MCM'
3079 ! include 'COMMON.MAP'
3080 ! include 'COMMON.HEADER'
3081 ! include 'COMMON.CSA'
3082 ! include 'COMMON.CHAIN'
3083 ! include 'COMMON.MUCA'
3084 ! include 'COMMON.MD'
3085 ! include 'COMMON.FFIELD'
3086 ! include 'COMMON.INTERACT'
3087 ! include 'COMMON.SETUP'
3088 !el integer :: KDIAG,ICORFL,IXDR
3089 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3090 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3091 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3092 ! character(len=80) :: ucase
3093 character(len=640) :: controlcard
3095 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3101 read (INP,'(a)') titel
3102 call card_concat(controlcard,.true.)
3103 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3104 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3105 call reada(controlcard,'SEED',seed,0.0D0)
3106 call random_init(seed)
3107 ! Set up the time limit (caution! The time must be input in minutes!)
3108 read_cart=index(controlcard,'READ_CART').gt.0
3109 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3110 call readi(controlcard,'SYM',symetr,1)
3111 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3112 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3113 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3114 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3115 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3116 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3117 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3118 call reada(controlcard,'DRMS',drms,0.1D0)
3119 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3120 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3121 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3122 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3123 write (iout,'(a,f10.1)')'DRMS = ',drms
3124 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3125 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3127 call readi(controlcard,'NZ_START',nz_start,0)
3128 call readi(controlcard,'NZ_END',nz_end,0)
3129 call readi(controlcard,'IZ_SC',iz_sc,0)
3130 timlim=60.0D0*timlim
3131 safety = 60.0d0*safety
3134 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3135 !C SHIELD keyword sets if the shielding effect of side-chains is used
3136 !C 0 denots no shielding is used all peptide are equally despite the
3137 !C solvent accesible area
3138 !C 1 the newly introduced function
3139 !C 2 reseved for further possible developement
3140 call readi(controlcard,'SHIELD',shield_mode,0)
3141 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3142 write(iout,*) "shield_mode",shield_mode
3143 !C Varibles set size of box
3144 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3145 protein=index(controlcard,"PROTEIN").gt.0
3146 ions=index(controlcard,"IONS").gt.0
3147 nucleic=index(controlcard,"NUCLEIC").gt.0
3148 write (iout,*) "with_theta_constr ",with_theta_constr
3149 AFMlog=(index(controlcard,'AFM'))
3150 selfguide=(index(controlcard,'SELFGUIDE'))
3151 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3152 call readi(controlcard,'GENCONSTR',genconstr,0)
3153 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3154 call reada(controlcard,'BOXY',boxysize,100.0d0)
3155 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3156 call readi(controlcard,'TUBEMOD',tubemode,0)
3157 write (iout,*) TUBEmode,"TUBEMODE"
3158 if (TUBEmode.gt.0) then
3159 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3160 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3161 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3162 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3163 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3164 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3165 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3166 buftubebot=bordtubebot+tubebufthick
3167 buftubetop=bordtubetop-tubebufthick
3170 ! CUTOFFF ON ELECTROSTATICS
3171 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3172 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3173 write(iout,*) "R_CUT_ELE=",r_cut_ele
3174 ! Lipidic parameters
3175 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3176 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3177 if (lipthick.gt.0.0d0) then
3178 bordliptop=(boxzsize+lipthick)/2.0
3179 bordlipbot=bordliptop-lipthick
3180 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3181 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3182 buflipbot=bordlipbot+lipbufthick
3183 bufliptop=bordliptop-lipbufthick
3184 if ((lipbufthick*2.0d0).gt.lipthick) &
3185 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3186 endif !lipthick.gt.0
3187 write(iout,*) "bordliptop=",bordliptop
3188 write(iout,*) "bordlipbot=",bordlipbot
3189 write(iout,*) "bufliptop=",bufliptop
3190 write(iout,*) "buflipbot=",buflipbot
3191 write (iout,*) "SHIELD MODE",shield_mode
3193 !C-------------------------
3194 minim=(index(controlcard,'MINIMIZE').gt.0)
3195 dccart=(index(controlcard,'CART').gt.0)
3196 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3197 overlapsc=.not.overlapsc
3198 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3199 searchsc=.not.searchsc
3200 sideadd=(index(controlcard,'SIDEADD').gt.0)
3201 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3202 outpdb=(index(controlcard,'PDBOUT').gt.0)
3203 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3204 pdbref=(index(controlcard,'PDBREF').gt.0)
3205 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3206 indpdb=index(controlcard,'PDBSTART')
3207 extconf=(index(controlcard,'EXTCONF').gt.0)
3208 call readi(controlcard,'IPRINT',iprint,0)
3209 call readi(controlcard,'MAXGEN',maxgen,10000)
3210 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3211 call readi(controlcard,"KDIAG",kdiag,0)
3212 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3213 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3214 write (iout,*) "RESCALE_MODE",rescale_mode
3215 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3216 if (index(controlcard,'REGULAR').gt.0.0D0) then
3217 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3221 if (index(controlcard,'CHECKGRAD').gt.0) then
3223 if (index(controlcard,'CART').gt.0) then
3225 elseif (index(controlcard,'CARINT').gt.0) then
3230 elseif (index(controlcard,'THREAD').gt.0) then
3232 call readi(controlcard,'THREAD',nthread,0)
3233 if (nthread.gt.0) then
3234 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3237 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3238 stop 'Error termination in Read_Control.'
3240 else if (index(controlcard,'MCMA').gt.0) then
3242 else if (index(controlcard,'MCEE').gt.0) then
3244 else if (index(controlcard,'MULTCONF').gt.0) then
3246 else if (index(controlcard,'MAP').gt.0) then
3248 call readi(controlcard,'MAP',nmap,0)
3249 else if (index(controlcard,'CSA').gt.0) then
3251 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3253 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3256 !fcm else if (index(controlcard,'MCMF').gt.0) then
3258 else if (index(controlcard,'SOFTREG').gt.0) then
3260 else if (index(controlcard,'CHECK_BOND').gt.0) then
3262 else if (index(controlcard,'TEST').gt.0) then
3264 else if (index(controlcard,'MD').gt.0) then
3266 else if (index(controlcard,'RE ').gt.0) then
3270 lmuca=index(controlcard,'MUCA').gt.0
3271 call readi(controlcard,'MUCADYN',mucadyn,0)
3272 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3273 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3275 write (iout,*) 'MUCADYN=',mucadyn
3276 write (iout,*) 'MUCASMOOTH=',muca_smooth
3279 iscode=index(controlcard,'ONE_LETTER')
3280 indphi=index(controlcard,'PHI')
3281 indback=index(controlcard,'BACK')
3282 iranconf=index(controlcard,'RAND_CONF')
3283 i2ndstr=index(controlcard,'USE_SEC_PRED')
3284 gradout=index(controlcard,'GRADOUT').gt.0
3285 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3286 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3287 if (me.eq.king .or. .not.out1file ) &
3288 write (iout,*) "DISTCHAINMAX",distchainmax
3290 if(me.eq.king.or..not.out1file) &
3291 write (iout,'(2a)') diagmeth(kdiag),&
3292 ' routine used to diagonalize matrices.'
3293 if (shield_mode.gt.0) then
3295 !C VSolvSphere the volume of solving sphere
3297 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3298 !C there will be no distinction between proline peptide group and normal peptide
3299 !C group in case of shielding parameters
3300 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3301 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3302 write (iout,*) VSolvSphere,VSolvSphere_div
3303 !C long axis of side chain
3305 long_r_sidechain(i)=vbldsc0(1,i)
3306 short_r_sidechain(i)=sigma0(i)
3307 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3313 end subroutine read_control
3314 !-----------------------------------------------------------------------------
3315 subroutine read_REMDpar
3317 ! Read REMD settings
3324 use control_data, only:out1file
3325 ! implicit real*8 (a-h,o-z)
3326 ! include 'DIMENSIONS'
3327 ! include 'COMMON.IOUNITS'
3328 ! include 'COMMON.TIME1'
3329 ! include 'COMMON.MD'
3332 !el include 'COMMON.LANGEVIN'
3334 !el include 'COMMON.LANGEVIN.lang0'
3336 ! include 'COMMON.INTERACT'
3337 ! include 'COMMON.NAMES'
3338 ! include 'COMMON.GEO'
3339 ! include 'COMMON.REMD'
3340 ! include 'COMMON.CONTROL'
3341 ! include 'COMMON.SETUP'
3342 ! character(len=80) :: ucase
3343 character(len=320) :: controlcard
3344 character(len=3200) :: controlcard1
3345 integer :: iremd_m_total
3348 ! real(kind=8) :: var,ene
3350 if(me.eq.king.or..not.out1file) &
3351 write (iout,*) "REMD setup"
3353 call card_concat(controlcard,.true.)
3354 call readi(controlcard,"NREP",nrep,3)
3355 call readi(controlcard,"NSTEX",nstex,1000)
3356 call reada(controlcard,"RETMIN",retmin,10.0d0)
3357 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3358 mremdsync=(index(controlcard,'SYNC').gt.0)
3359 call readi(controlcard,"NSYN",i_sync_step,100)
3360 restart1file=(index(controlcard,'REST1FILE').gt.0)
3361 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3362 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3363 if(max_cache_traj_use.gt.max_cache_traj) &
3364 max_cache_traj_use=max_cache_traj
3365 if(me.eq.king.or..not.out1file) then
3366 !d if (traj1file) then
3367 !rc caching is in testing - NTWX is not ignored
3368 !d write (iout,*) "NTWX value is ignored"
3369 !d write (iout,*) " trajectory is stored to one file by master"
3370 !d write (iout,*) " before exchange at NSTEX intervals"
3372 write (iout,*) "NREP= ",nrep
3373 write (iout,*) "NSTEX= ",nstex
3374 write (iout,*) "SYNC= ",mremdsync
3375 write (iout,*) "NSYN= ",i_sync_step
3376 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3379 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3380 if (index(controlcard,'TLIST').gt.0) then
3382 call card_concat(controlcard1,.true.)
3383 read(controlcard1,*) (remd_t(i),i=1,nrep)
3384 if(me.eq.king.or..not.out1file) &
3385 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3388 if (index(controlcard,'MLIST').gt.0) then
3390 call card_concat(controlcard1,.true.)
3391 read(controlcard1,*) (remd_m(i),i=1,nrep)
3392 if(me.eq.king.or..not.out1file) then
3393 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3396 iremd_m_total=iremd_m_total+remd_m(i)
3398 write (iout,*) 'Total number of replicas ',iremd_m_total
3401 if(me.eq.king.or..not.out1file) &
3402 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3404 end subroutine read_REMDpar
3405 !-----------------------------------------------------------------------------
3406 subroutine read_MDpar
3410 use control_data, only: r_cut,rlamb,out1file
3412 use geometry_data, only: pi
3414 ! implicit real*8 (a-h,o-z)
3415 ! include 'DIMENSIONS'
3416 ! include 'COMMON.IOUNITS'
3417 ! include 'COMMON.TIME1'
3418 ! include 'COMMON.MD'
3421 !el include 'COMMON.LANGEVIN'
3423 !el include 'COMMON.LANGEVIN.lang0'
3425 ! include 'COMMON.INTERACT'
3426 ! include 'COMMON.NAMES'
3427 ! include 'COMMON.GEO'
3428 ! include 'COMMON.SETUP'
3429 ! include 'COMMON.CONTROL'
3430 ! include 'COMMON.SPLITELE'
3431 ! character(len=80) :: ucase
3432 character(len=320) :: controlcard
3437 call card_concat(controlcard,.true.)
3438 call readi(controlcard,"NSTEP",n_timestep,1000000)
3439 call readi(controlcard,"NTWE",ntwe,100)
3440 call readi(controlcard,"NTWX",ntwx,1000)
3441 call reada(controlcard,"DT",d_time,1.0d-1)
3442 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3443 call reada(controlcard,"DAMAX",damax,1.0d1)
3444 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3445 call readi(controlcard,"LANG",lang,0)
3446 RESPA = index(controlcard,"RESPA") .gt. 0
3447 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3448 ntime_split0=ntime_split
3449 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3450 ntime_split0=ntime_split
3451 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3452 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3453 rest = index(controlcard,"REST").gt.0
3454 tbf = index(controlcard,"TBF").gt.0
3455 usampl = index(controlcard,"USAMPL").gt.0
3456 mdpdb = index(controlcard,"MDPDB").gt.0
3457 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3458 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3459 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3460 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3461 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3462 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3463 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3464 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3465 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3466 large = index(controlcard,"LARGE").gt.0
3467 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3468 rattle = index(controlcard,"RATTLE").gt.0
3469 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3475 if(me.eq.king.or..not.out1file) then
3477 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3479 write (iout,'(a)') "The units are:"
3480 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3481 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3482 " acceleration: angstrom/(48.9 fs)**2"
3483 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3485 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3486 write (iout,'(a60,f10.5,a)') &
3487 "Initial time step of numerical integration:",d_time,&
3489 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3491 write (iout,'(2a,i4,a)') &
3492 "A-MTS algorithm used; initial time step for fast-varying",&
3493 " short-range forces split into",ntime_split," steps."
3494 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3495 r_cut," lambda",rlamb
3497 write (iout,'(2a,f10.5)') &
3498 "Maximum acceleration threshold to reduce the time step",&
3499 "/increase split number:",damax
3500 write (iout,'(2a,f10.5)') &
3501 "Maximum predicted energy drift to reduce the timestep",&
3502 "/increase split number:",edriftmax
3503 write (iout,'(a60,f10.5)') &
3504 "Maximum velocity threshold to reduce velocities:",dvmax
3505 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3506 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3507 if (rattle) write (iout,'(a60)') &
3508 "Rattle algorithm used to constrain the virtual bonds"
3512 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3513 call reada(controlcard,"RWAT",rwat,1.4d0)
3514 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3515 surfarea=index(controlcard,"SURFAREA").gt.0
3516 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3517 if(me.eq.king.or..not.out1file)then
3518 write (iout,'(/a,$)') "Langevin dynamics calculation"
3520 write (iout,'(a/)') &
3521 " with direct integration of Langevin equations"
3522 else if (lang.eq.2) then
3523 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3524 else if (lang.eq.3) then
3525 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3526 else if (lang.eq.4) then
3527 write (iout,'(a/)') " in overdamped mode"
3529 write (iout,'(//a,i5)') &
3530 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3533 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3534 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3535 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3536 write (iout,'(a60,f10.5)') &
3537 "Scaling factor of the friction forces:",scal_fric
3538 if (surfarea) write (iout,'(2a,i10,a)') &
3539 "Friction coefficients will be scaled by solvent-accessible",&
3540 " surface area every",reset_fricmat," steps."
3542 ! Calculate friction coefficients and bounds of stochastic forces
3543 eta=6*pi*cPoise*etawat
3544 if(me.eq.king.or..not.out1file) &
3545 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3547 gamp=scal_fric*(pstok+rwat)*eta
3548 stdfp=dsqrt(2*Rb*t_bath/d_time)
3549 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3551 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3552 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3554 if(me.eq.king.or..not.out1file)then
3555 write (iout,'(/2a/)') &
3556 "Radii of site types and friction coefficients and std's of",&
3557 " stochastic forces of fully exposed sites"
3558 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3560 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i),&
3561 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3565 if(me.eq.king.or..not.out1file)then
3566 write (iout,'(a)') "Berendsen bath calculation"
3567 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3568 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3570 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3571 count_reset_moment," steps"
3573 write (iout,'(a,i10,a)') &
3574 "Velocities will be reset at random every",count_reset_vel,&
3578 if(me.eq.king.or..not.out1file) &
3579 write (iout,'(a31)') "Microcanonical mode calculation"
3581 if(me.eq.king.or..not.out1file)then
3582 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3584 write(iout,*) "MD running with constraints."
3585 write(iout,*) "Equilibration time ", eq_time, " mtus."
3586 write(iout,*) "Constraining ", nfrag," fragments."
3587 write(iout,*) "Length of each fragment, weight and q0:"
3589 write (iout,*) "Set of restraints #",iset
3591 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3592 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3594 write(iout,*) "constraints between ", npair, "fragments."
3595 write(iout,*) "constraint pairs, weights and q0:"
3597 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3598 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3600 write(iout,*) "angle constraints within ", nfrag_back,&
3601 "backbone fragments."
3602 write(iout,*) "fragment, weights:"
3604 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3605 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3606 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3609 iset=mod(kolor,nset)+1
3612 if(me.eq.king.or..not.out1file) &
3613 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3615 end subroutine read_MDpar
3616 !-----------------------------------------------------------------------------
3620 ! implicit real*8 (a-h,o-z)
3621 ! include 'DIMENSIONS'
3622 ! include 'COMMON.MAP'
3623 ! include 'COMMON.IOUNITS'
3624 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3625 character(len=80) :: mapcard !,ucase
3628 ! real(kind=8) :: var,ene
3631 read (inp,'(a)') mapcard
3632 mapcard=ucase(mapcard)
3633 if (index(mapcard,'PHI').gt.0) then
3635 else if (index(mapcard,'THE').gt.0) then
3637 else if (index(mapcard,'ALP').gt.0) then
3639 else if (index(mapcard,'OME').gt.0) then
3642 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3643 stop 'Error - illegal variable spec in MAP card.'
3645 call readi (mapcard,'RES1',res1(imap),0)
3646 call readi (mapcard,'RES2',res2(imap),0)
3647 if (res1(imap).eq.0) then
3648 res1(imap)=res2(imap)
3649 else if (res2(imap).eq.0) then
3650 res2(imap)=res1(imap)
3652 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3653 write (iout,'(a)') &
3654 'Error - illegal definition of variable group in MAP.'
3655 stop 'Error - illegal definition of variable group in MAP.'
3657 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3658 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3659 call readi(mapcard,'NSTEP',nstep(imap),0)
3660 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3661 write (iout,'(a)') &
3662 'Illegal boundary and/or step size specification in MAP.'
3663 stop 'Illegal boundary and/or step size specification in MAP.'
3667 end subroutine map_read
3668 !-----------------------------------------------------------------------------
3671 use control_data, only: vdisulf
3673 ! implicit real*8 (a-h,o-z)
3674 ! include 'DIMENSIONS'
3675 ! include 'COMMON.IOUNITS'
3676 ! include 'COMMON.GEO'
3677 ! include 'COMMON.CSA'
3678 ! include 'COMMON.BANK'
3679 ! include 'COMMON.CONTROL'
3680 ! character(len=80) :: ucase
3681 character(len=620) :: mcmcard
3683 ! integer :: ntf,ik,iw_pdb
3684 ! real(kind=8) :: var,ene
3686 call card_concat(mcmcard,.true.)
3688 call readi(mcmcard,'NCONF',nconf,50)
3689 call readi(mcmcard,'NADD',nadd,0)
3690 call readi(mcmcard,'JSTART',jstart,1)
3691 call readi(mcmcard,'JEND',jend,1)
3692 call readi(mcmcard,'NSTMAX',nstmax,500000)
3693 call readi(mcmcard,'N0',n0,1)
3694 call readi(mcmcard,'N1',n1,6)
3695 call readi(mcmcard,'N2',n2,4)
3696 call readi(mcmcard,'N3',n3,0)
3697 call readi(mcmcard,'N4',n4,0)
3698 call readi(mcmcard,'N5',n5,0)
3699 call readi(mcmcard,'N6',n6,10)
3700 call readi(mcmcard,'N7',n7,0)
3701 call readi(mcmcard,'N8',n8,0)
3702 call readi(mcmcard,'N9',n9,0)
3703 call readi(mcmcard,'N14',n14,0)
3704 call readi(mcmcard,'N15',n15,0)
3705 call readi(mcmcard,'N16',n16,0)
3706 call readi(mcmcard,'N17',n17,0)
3707 call readi(mcmcard,'N18',n18,0)
3709 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3711 call readi(mcmcard,'NDIFF',ndiff,2)
3712 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3713 call readi(mcmcard,'IS1',is1,1)
3714 call readi(mcmcard,'IS2',is2,8)
3715 call readi(mcmcard,'NRAN0',nran0,4)
3716 call readi(mcmcard,'NRAN1',nran1,2)
3717 call readi(mcmcard,'IRR',irr,1)
3718 call readi(mcmcard,'NSEED',nseed,20)
3719 call readi(mcmcard,'NTOTAL',ntotal,10000)
3720 call reada(mcmcard,'CUT1',cut1,2.0d0)
3721 call reada(mcmcard,'CUT2',cut2,5.0d0)
3722 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3723 call readi(mcmcard,'ICMAX',icmax,3)
3724 call readi(mcmcard,'IRESTART',irestart,0)
3725 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3728 call reada(mcmcard,'DELE',dele,20.0d0)
3729 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3730 call readi(mcmcard,'IREF',iref,0)
3731 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3732 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3733 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3734 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3735 write (iout,*) "NCONF_IN",nconf_in
3737 end subroutine csaread
3738 !-----------------------------------------------------------------------------
3742 use control_data, only: MaxMoveType
3745 ! implicit real*8 (a-h,o-z)
3746 ! include 'DIMENSIONS'
3747 ! include 'COMMON.MCM'
3748 ! include 'COMMON.MCE'
3749 ! include 'COMMON.IOUNITS'
3750 ! character(len=80) :: ucase
3751 character(len=320) :: mcmcard
3754 ! real(kind=8) :: var,ene
3756 call card_concat(mcmcard,.true.)
3757 call readi(mcmcard,'MAXACC',maxacc,100)
3758 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3759 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3760 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3761 call readi(mcmcard,'MAXREPM',maxrepm,200)
3762 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3763 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3764 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3765 call reada(mcmcard,'E_UP',e_up,5.0D0)
3766 call reada(mcmcard,'DELTE',delte,0.1D0)
3767 call readi(mcmcard,'NSWEEP',nsweep,5)
3768 call readi(mcmcard,'NSTEPH',nsteph,0)
3769 call readi(mcmcard,'NSTEPC',nstepc,0)
3770 call reada(mcmcard,'TMIN',tmin,298.0D0)
3771 call reada(mcmcard,'TMAX',tmax,298.0D0)
3772 call readi(mcmcard,'NWINDOW',nwindow,0)
3773 call readi(mcmcard,'PRINT_MC',print_mc,0)
3774 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3775 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3776 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3777 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3778 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3779 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3780 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3781 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3782 if (nwindow.gt.0) then
3783 allocate(winstart(nwindow)) !!el (maxres)
3784 allocate(winend(nwindow)) !!el
3785 allocate(winlen(nwindow)) !!el
3786 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3788 winlen(i)=winend(i)-winstart(i)+1
3791 if (tmax.lt.tmin) tmax=tmin
3792 if (tmax.eq.tmin) then
3796 if (nstepc.gt.0 .and. nsteph.gt.0) then
3797 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3798 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3800 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3801 ! Probabilities of different move types
3802 sumpro_type(0)=0.0D0
3803 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3804 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3805 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3806 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3807 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3808 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3809 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3811 print *,'i',i,' sumprotype',sumpro_type(i)
3812 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3813 print *,'i',i,' sumprotype',sumpro_type(i)
3816 end subroutine mcmread
3817 !-----------------------------------------------------------------------------
3818 subroutine read_minim
3821 ! implicit real*8 (a-h,o-z)
3822 ! include 'DIMENSIONS'
3823 ! include 'COMMON.MINIM'
3824 ! include 'COMMON.IOUNITS'
3825 ! character(len=80) :: ucase
3826 character(len=320) :: minimcard
3828 ! integer :: ntf,ik,iw_pdb
3829 ! real(kind=8) :: var,ene
3831 call card_concat(minimcard,.true.)
3832 call readi(minimcard,'MAXMIN',maxmin,2000)
3833 call readi(minimcard,'MAXFUN',maxfun,5000)
3834 call readi(minimcard,'MINMIN',minmin,maxmin)
3835 call readi(minimcard,'MINFUN',minfun,maxmin)
3836 call reada(minimcard,'TOLF',tolf,1.0D-2)
3837 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3838 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3839 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3840 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3841 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3842 'Options in energy minimization:'
3843 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3844 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3845 'MinMin:',MinMin,' MinFun:',MinFun,&
3846 ' TolF:',TolF,' RTolF:',RTolF
3848 end subroutine read_minim
3849 !-----------------------------------------------------------------------------
3850 subroutine openunits
3852 use MD_data, only: usampl
3855 use control_data, only:out1file
3856 use control, only: getenv_loc
3857 ! implicit real*8 (a-h,o-z)
3858 ! include 'DIMENSIONS'
3861 character(len=16) :: form,nodename
3862 integer :: nodelen,ierror,npos
3864 ! include 'COMMON.SETUP'
3865 ! include 'COMMON.IOUNITS'
3866 ! include 'COMMON.MD'
3867 ! include 'COMMON.CONTROL'
3868 integer :: lenpre,lenpot,lentmp !,ilen
3870 character(len=3) :: out1file_text !,ucase
3871 character(len=3) :: ll
3874 ! integer :: ntf,ik,iw_pdb
3875 ! real(kind=8) :: var,ene
3877 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3878 call getenv_loc("PREFIX",prefix)
3880 call getenv_loc("POT",pot)
3881 call getenv_loc("DIRTMP",tmpdir)
3882 call getenv_loc("CURDIR",curdir)
3883 call getenv_loc("OUT1FILE",out1file_text)
3884 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3885 out1file_text=ucase(out1file_text)
3886 if (out1file_text(1:1).eq."Y") then
3889 out1file=fg_rank.gt.0
3894 if (lentmp.gt.0) then
3895 write (*,'(80(1h!))')
3896 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3897 write (*,'(80(1h!))')
3898 write (*,*)"All output files will be on node /tmp directory."
3900 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3901 if (me.eq.king) then
3902 write (*,*) "The master node is ",nodename
3903 else if (fg_rank.eq.0) then
3904 write (*,*) "I am the CG slave node ",nodename
3906 write (*,*) "I am the FG slave node ",nodename
3909 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3910 lenpre = lentmp+lenpre+1
3912 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3913 ! Get the names and open the input files
3914 #if defined(WINIFL) || defined(WINPGI)
3915 open(1,file=pref_orig(:ilen(pref_orig))// &
3916 '.inp',status='old',readonly,shared)
3917 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3918 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3919 ! Get parameter filenames and open the parameter files.
3920 call getenv_loc('BONDPAR',bondname)
3921 open (ibond,file=bondname,status='old',readonly,shared)
3922 call getenv_loc('THETPAR',thetname)
3923 open (ithep,file=thetname,status='old',readonly,shared)
3924 call getenv_loc('ROTPAR',rotname)
3925 open (irotam,file=rotname,status='old',readonly,shared)
3926 call getenv_loc('TORPAR',torname)
3927 open (itorp,file=torname,status='old',readonly,shared)
3928 call getenv_loc('TORDPAR',tordname)
3929 open (itordp,file=tordname,status='old',readonly,shared)
3930 call getenv_loc('FOURIER',fouriername)
3931 open (ifourier,file=fouriername,status='old',readonly,shared)
3932 call getenv_loc('ELEPAR',elename)
3933 open (ielep,file=elename,status='old',readonly,shared)
3934 call getenv_loc('SIDEPAR',sidename)
3935 open (isidep,file=sidename,status='old',readonly,shared)
3936 #elif (defined CRAY) || (defined AIX)
3937 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3939 ! print *,"Processor",myrank," opened file 1"
3940 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3941 ! print *,"Processor",myrank," opened file 9"
3942 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3943 ! Get parameter filenames and open the parameter files.
3944 call getenv_loc('BONDPAR',bondname)
3945 open (ibond,file=bondname,status='old',action='read')
3946 ! print *,"Processor",myrank," opened file IBOND"
3947 call getenv_loc('THETPAR',thetname)
3948 open (ithep,file=thetname,status='old',action='read')
3949 ! print *,"Processor",myrank," opened file ITHEP"
3950 call getenv_loc('ROTPAR',rotname)
3951 open (irotam,file=rotname,status='old',action='read')
3952 ! print *,"Processor",myrank," opened file IROTAM"
3953 call getenv_loc('TORPAR',torname)
3954 open (itorp,file=torname,status='old',action='read')
3955 ! print *,"Processor",myrank," opened file ITORP"
3956 call getenv_loc('TORDPAR',tordname)
3957 open (itordp,file=tordname,status='old',action='read')
3958 ! print *,"Processor",myrank," opened file ITORDP"
3959 call getenv_loc('SCCORPAR',sccorname)
3960 open (isccor,file=sccorname,status='old',action='read')
3961 ! print *,"Processor",myrank," opened file ISCCOR"
3962 call getenv_loc('FOURIER',fouriername)
3963 open (ifourier,file=fouriername,status='old',action='read')
3964 ! print *,"Processor",myrank," opened file IFOURIER"
3965 call getenv_loc('ELEPAR',elename)
3966 open (ielep,file=elename,status='old',action='read')
3967 ! print *,"Processor",myrank," opened file IELEP"
3968 call getenv_loc('SIDEPAR',sidename)
3969 open (isidep,file=sidename,status='old',action='read')
3970 call getenv_loc('LIPTRANPAR',liptranname)
3971 open (iliptranpar,file=liptranname,status='old',action='read')
3972 call getenv_loc('TUBEPAR',tubename)
3973 open (itube,file=tubename,status='old',action='read')
3975 ! print *,"Processor",myrank," opened file ISIDEP"
3976 ! print *,"Processor",myrank," opened parameter files"
3978 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3979 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3980 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3981 ! Get parameter filenames and open the parameter files.
3982 call getenv_loc('BONDPAR',bondname)
3983 open (ibond,file=bondname,status='old')
3984 call getenv_loc('THETPAR',thetname)
3985 open (ithep,file=thetname,status='old')
3986 call getenv_loc('ROTPAR',rotname)
3987 open (irotam,file=rotname,status='old')
3988 call getenv_loc('TORPAR',torname)
3989 open (itorp,file=torname,status='old')
3990 call getenv_loc('TORDPAR',tordname)
3991 open (itordp,file=tordname,status='old')
3992 call getenv_loc('SCCORPAR',sccorname)
3993 open (isccor,file=sccorname,status='old')
3994 call getenv_loc('FOURIER',fouriername)
3995 open (ifourier,file=fouriername,status='old')
3996 call getenv_loc('ELEPAR',elename)
3997 open (ielep,file=elename,status='old')
3998 call getenv_loc('SIDEPAR',sidename)
3999 open (isidep,file=sidename,status='old')
4000 call getenv_loc('LIPTRANPAR',liptranname)
4001 open (iliptranpar,file=liptranname,status='old')
4002 call getenv_loc('TUBEPAR',tubename)
4003 open (itube,file=tubename,status='old')
4005 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4007 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4008 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4009 ! Get parameter filenames and open the parameter files.
4010 call getenv_loc('BONDPAR',bondname)
4011 open (ibond,file=bondname,status='old',action='read')
4012 call getenv_loc('THETPAR',thetname)
4013 open (ithep,file=thetname,status='old',action='read')
4014 call getenv_loc('ROTPAR',rotname)
4015 open (irotam,file=rotname,status='old',action='read')
4016 call getenv_loc('TORPAR',torname)
4017 open (itorp,file=torname,status='old',action='read')
4018 call getenv_loc('TORDPAR',tordname)
4019 open (itordp,file=tordname,status='old',action='read')
4020 call getenv_loc('SCCORPAR',sccorname)
4021 open (isccor,file=sccorname,status='old',action='read')
4023 call getenv_loc('THETPARPDB',thetname_pdb)
4024 print *,"thetname_pdb ",thetname_pdb
4025 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4026 print *,ithep_pdb," opened"
4028 call getenv_loc('FOURIER',fouriername)
4029 open (ifourier,file=fouriername,status='old',readonly)
4030 call getenv_loc('ELEPAR',elename)
4031 open (ielep,file=elename,status='old',readonly)
4032 call getenv_loc('SIDEPAR',sidename)
4033 open (isidep,file=sidename,status='old',readonly)
4034 call getenv_loc('LIPTRANPAR',liptranname)
4035 open (iliptranpar,file=liptranname,status='old',action='read')
4036 call getenv_loc('TUBEPAR',tubename)
4037 open (itube,file=tubename,status='old',action='read')
4040 call getenv_loc('ROTPARPDB',rotname_pdb)
4041 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4046 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4047 ! Use -DOLDSCP to use hard-coded constants instead.
4049 call getenv_loc('SCPPAR',scpname)
4050 #if defined(WINIFL) || defined(WINPGI)
4051 open (iscpp,file=scpname,status='old',readonly,shared)
4052 #elif (defined CRAY) || (defined AIX)
4053 open (iscpp,file=scpname,status='old',action='read')
4055 open (iscpp,file=scpname,status='old')
4057 open (iscpp,file=scpname,status='old',action='read')
4060 call getenv_loc('PATTERN',patname)
4061 #if defined(WINIFL) || defined(WINPGI)
4062 open (icbase,file=patname,status='old',readonly,shared)
4063 #elif (defined CRAY) || (defined AIX)
4064 open (icbase,file=patname,status='old',action='read')
4066 open (icbase,file=patname,status='old')
4068 open (icbase,file=patname,status='old',action='read')
4071 ! Open output file only for CG processes
4072 ! print *,"Processor",myrank," fg_rank",fg_rank
4073 if (fg_rank.eq.0) then
4075 if (nodes.eq.1) then
4078 npos = dlog10(dfloat(nodes-1))+1
4080 if (npos.lt.3) npos=3
4081 write (liczba,'(i1)') npos
4082 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4084 write (liczba,form) me
4085 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4086 liczba(:ilen(liczba))
4087 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4089 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4091 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4092 liczba(:ilen(liczba))//'.mol2'
4093 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4094 liczba(:ilen(liczba))//'.stat'
4096 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4097 //liczba(:ilen(liczba))//'.stat')
4098 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4101 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4102 liczba(:ilen(liczba))//'.const'
4107 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4108 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4109 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4110 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4111 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4113 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4115 rest2name=prefix(:ilen(prefix))//'.rst'
4117 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4120 #if defined(AIX) || defined(PGI)
4121 if (me.eq.king .or. .not. out1file) &
4122 open(iout,file=outname,status='unknown')
4124 if (fg_rank.gt.0) then
4125 write (liczba,'(i3.3)') myrank/nfgtasks
4126 write (ll,'(bz,i3.3)') fg_rank
4127 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4132 open(igeom,file=intname,status='unknown',position='append')
4133 open(ipdb,file=pdbname,status='unknown')
4134 open(imol2,file=mol2name,status='unknown')
4135 open(istat,file=statname,status='unknown',position='append')
4137 !1out open(iout,file=outname,status='unknown')
4140 if (me.eq.king .or. .not.out1file) &
4141 open(iout,file=outname,status='unknown')
4143 if (fg_rank.gt.0) then
4144 write (liczba,'(i3.3)') myrank/nfgtasks
4145 write (ll,'(bz,i3.3)') fg_rank
4146 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4151 open(igeom,file=intname,status='unknown',access='append')
4152 open(ipdb,file=pdbname,status='unknown')
4153 open(imol2,file=mol2name,status='unknown')
4154 open(istat,file=statname,status='unknown',access='append')
4156 !1out open(iout,file=outname,status='unknown')
4159 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4160 csa_seed=prefix(:lenpre)//'.CSA.seed'
4161 csa_history=prefix(:lenpre)//'.CSA.history'
4162 csa_bank=prefix(:lenpre)//'.CSA.bank'
4163 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4164 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4165 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4166 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4167 csa_int=prefix(:lenpre)//'.int'
4168 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4169 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4170 csa_in=prefix(:lenpre)//'.CSA.in'
4171 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4174 write (iout,'(80(1h-))')
4175 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4176 write (iout,'(80(1h-))')
4177 write (iout,*) "Input file : ",&
4178 pref_orig(:ilen(pref_orig))//'.inp'
4179 write (iout,*) "Output file : ",&
4180 outname(:ilen(outname))
4182 write (iout,*) "Sidechain potential file : ",&
4183 sidename(:ilen(sidename))
4185 write (iout,*) "SCp potential file : ",&
4186 scpname(:ilen(scpname))
4188 write (iout,*) "Electrostatic potential file : ",&
4189 elename(:ilen(elename))
4190 write (iout,*) "Cumulant coefficient file : ",&
4191 fouriername(:ilen(fouriername))
4192 write (iout,*) "Torsional parameter file : ",&
4193 torname(:ilen(torname))
4194 write (iout,*) "Double torsional parameter file : ",&
4195 tordname(:ilen(tordname))
4196 write (iout,*) "SCCOR parameter file : ",&
4197 sccorname(:ilen(sccorname))
4198 write (iout,*) "Bond & inertia constant file : ",&
4199 bondname(:ilen(bondname))
4200 write (iout,*) "Bending parameter file : ",&
4201 thetname(:ilen(thetname))
4202 write (iout,*) "Rotamer parameter file : ",&
4203 rotname(:ilen(rotname))
4206 write (iout,*) "Thetpdb parameter file : ",&
4207 thetname_pdb(:ilen(thetname_pdb))
4210 write (iout,*) "Threading database : ",&
4211 patname(:ilen(patname))
4213 write (iout,*)" DIRTMP : ",&
4215 write (iout,'(80(1h-))')
4218 end subroutine openunits
4219 !-----------------------------------------------------------------------------
4222 use geometry_data, only: nres,dc
4224 ! implicit real*8 (a-h,o-z)
4225 ! include 'DIMENSIONS'
4226 ! include 'COMMON.CHAIN'
4227 ! include 'COMMON.IOUNITS'
4228 ! include 'COMMON.MD'
4231 ! real(kind=8) :: var,ene
4233 open(irest2,file=rest2name,status='unknown')
4234 read(irest2,*) totT,EK,potE,totE,t_bath
4237 ! AL 4/17/17: Now reading d_t(0,:) too
4239 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4242 ! AL 4/17/17: Now reading d_c(0,:) too
4244 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4247 read (irest2,*) iset
4251 end subroutine readrst
4252 !-----------------------------------------------------------------------------
4253 subroutine read_fragments
4257 use control_data, only:out1file
4260 ! implicit real*8 (a-h,o-z)
4261 ! include 'DIMENSIONS'
4265 ! include 'COMMON.SETUP'
4266 ! include 'COMMON.CHAIN'
4267 ! include 'COMMON.IOUNITS'
4268 ! include 'COMMON.MD'
4269 ! include 'COMMON.CONTROL'
4272 ! real(kind=8) :: var,ene
4274 read(inp,*) nset,nfrag,npair,nfrag_back
4276 !el from module energy
4277 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4278 if(.not.allocated(wfrag_back)) then
4279 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4280 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4282 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4283 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4285 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4286 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4289 if(me.eq.king.or..not.out1file) &
4290 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4291 " nfrag_back",nfrag_back
4293 read(inp,*) mset(iset)
4295 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4297 if(me.eq.king.or..not.out1file) &
4298 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4299 ifrag(2,i,iset), qinfrag(i,iset)
4302 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4304 if(me.eq.king.or..not.out1file) &
4305 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4306 ipair(2,i,iset), qinpair(i,iset)
4309 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4310 wfrag_back(3,i,iset),&
4311 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4312 if(me.eq.king.or..not.out1file) &
4313 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4314 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4318 end subroutine read_fragments
4319 !-----------------------------------------------------------------------------
4321 !-----------------------------------------------------------------------------
4325 ! implicit real*8 (a-h,o-z)
4326 ! include 'DIMENSIONS'
4327 ! include 'COMMON.CSA'
4328 ! include 'COMMON.BANK'
4329 ! include 'COMMON.IOUNITS'
4331 ! integer :: ntf,ik,iw_pdb
4332 ! real(kind=8) :: var,ene
4334 open(icsa_in,file=csa_in,status="old",err=100)
4335 read(icsa_in,*) nconf
4336 read(icsa_in,*) jstart,jend
4337 read(icsa_in,*) nstmax
4338 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4339 read(icsa_in,*) nran0,nran1,irr
4340 read(icsa_in,*) nseed
4341 read(icsa_in,*) ntotal,cut1,cut2
4342 read(icsa_in,*) estop
4343 read(icsa_in,*) icmax,irestart
4344 read(icsa_in,*) ntbankm,dele,difcut
4345 read(icsa_in,*) iref,rmscut,pnccut
4346 read(icsa_in,*) ndiff
4353 end subroutine csa_read
4354 !-----------------------------------------------------------------------------
4355 subroutine initial_write
4358 ! implicit real*8 (a-h,o-z)
4359 ! include 'DIMENSIONS'
4360 ! include 'COMMON.CSA'
4361 ! include 'COMMON.BANK'
4362 ! include 'COMMON.IOUNITS'
4364 ! integer :: ntf,ik,iw_pdb
4365 ! real(kind=8) :: var,ene
4367 open(icsa_seed,file=csa_seed,status="unknown")
4368 write(icsa_seed,*) "seed"
4370 #if defined(AIX) || defined(PGI)
4371 open(icsa_history,file=csa_history,status="unknown",&
4374 open(icsa_history,file=csa_history,status="unknown",&
4377 write(icsa_history,*) nconf
4378 write(icsa_history,*) jstart,jend
4379 write(icsa_history,*) nstmax
4380 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4381 write(icsa_history,*) nran0,nran1,irr
4382 write(icsa_history,*) nseed
4383 write(icsa_history,*) ntotal,cut1,cut2
4384 write(icsa_history,*) estop
4385 write(icsa_history,*) icmax,irestart
4386 write(icsa_history,*) ntbankm,dele,difcut
4387 write(icsa_history,*) iref,rmscut,pnccut
4388 write(icsa_history,*) ndiff
4390 write(icsa_history,*)
4393 open(icsa_bank1,file=csa_bank1,status="unknown")
4394 write(icsa_bank1,*) 0
4398 end subroutine initial_write
4399 !-----------------------------------------------------------------------------
4400 subroutine restart_write
4403 ! implicit real*8 (a-h,o-z)
4404 ! include 'DIMENSIONS'
4405 ! include 'COMMON.IOUNITS'
4406 ! include 'COMMON.CSA'
4407 ! include 'COMMON.BANK'
4409 ! integer :: ntf,ik,iw_pdb
4410 ! real(kind=8) :: var,ene
4412 #if defined(AIX) || defined(PGI)
4413 open(icsa_history,file=csa_history,position="append")
4415 open(icsa_history,file=csa_history,access="append")
4417 write(icsa_history,*)
4418 write(icsa_history,*) "This is restart"
4419 write(icsa_history,*)
4420 write(icsa_history,*) nconf
4421 write(icsa_history,*) jstart,jend
4422 write(icsa_history,*) nstmax
4423 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4424 write(icsa_history,*) nran0,nran1,irr
4425 write(icsa_history,*) nseed
4426 write(icsa_history,*) ntotal,cut1,cut2
4427 write(icsa_history,*) estop
4428 write(icsa_history,*) icmax,irestart
4429 write(icsa_history,*) ntbankm,dele,difcut
4430 write(icsa_history,*) iref,rmscut,pnccut
4431 write(icsa_history,*) ndiff
4432 write(icsa_history,*)
4433 write(icsa_history,*) "irestart is: ", irestart
4435 write(icsa_history,*)
4439 end subroutine restart_write
4440 !-----------------------------------------------------------------------------
4442 !-----------------------------------------------------------------------------
4443 subroutine write_pdb(npdb,titelloc,ee)
4445 ! implicit real*8 (a-h,o-z)
4446 ! include 'DIMENSIONS'
4447 ! include 'COMMON.IOUNITS'
4448 character(len=50) :: titelloc1
4449 character*(*) :: titelloc
4450 character(len=3) :: zahl
4451 character(len=5) :: liczba5
4453 integer :: npdb !,ilen
4457 ! real(kind=8) :: var,ene
4461 if (npdb.lt.1000) then
4462 call numstr(npdb,zahl)
4463 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4465 if (npdb.lt.10000) then
4466 write(liczba5,'(i1,i4)') 0,npdb
4468 write(liczba5,'(i5)') npdb
4470 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4472 call pdbout(ee,titelloc1,ipdb)
4475 end subroutine write_pdb
4476 !-----------------------------------------------------------------------------
4478 !-----------------------------------------------------------------------------
4479 subroutine write_thread_summary
4480 ! Thread the sequence through a database of known structures
4481 use control_data, only: refstr
4483 use energy_data, only: n_ene_comp
4485 ! implicit real*8 (a-h,o-z)
4486 ! include 'DIMENSIONS'
4488 use MPI_data !include 'COMMON.INFO'
4491 ! include 'COMMON.CONTROL'
4492 ! include 'COMMON.CHAIN'
4493 ! include 'COMMON.DBASE'
4494 ! include 'COMMON.INTERACT'
4495 ! include 'COMMON.VAR'
4496 ! include 'COMMON.THREAD'
4497 ! include 'COMMON.FFIELD'
4498 ! include 'COMMON.SBRIDGE'
4499 ! include 'COMMON.HEADER'
4500 ! include 'COMMON.NAMES'
4501 ! include 'COMMON.IOUNITS'
4502 ! include 'COMMON.TIME1'
4504 integer,dimension(maxthread) :: ip
4505 real(kind=8),dimension(0:n_ene) :: energia
4507 integer :: i,j,ii,jj,ipj,ik,kk,ist
4508 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4510 write (iout,'(30x,a/)') &
4511 ' *********** Summary threading statistics ************'
4512 write (iout,'(a)') 'Initial energies:'
4513 write (iout,'(a4,2x,a12,14a14,3a8)') &
4514 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4515 'RMSnat','NatCONT','NNCONT','RMS'
4516 ! Energy sort patterns
4521 enet=ener(n_ene-1,ip(i))
4524 if (ener(n_ene-1,ip(j)).lt.enet) then
4526 enet=ener(n_ene-1,ip(j))
4538 ist=nres_base(2,ii)+ipatt(2,i)
4540 energia(i)=ener0(kk,i)
4542 etot=ener0(n_ene_comp+1,i)
4543 rmsnat=ener0(n_ene_comp+2,i)
4544 rms=ener0(n_ene_comp+3,i)
4545 frac=ener0(n_ene_comp+4,i)
4546 frac_nn=ener0(n_ene_comp+5,i)
4549 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4550 i,str_nam(ii),ist+1,&
4551 (energia(print_order(kk)),kk=1,nprint_ene),&
4552 etot,rmsnat,frac,frac_nn,rms
4554 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4555 i,str_nam(ii),ist+1,&
4556 (energia(print_order(kk)),kk=1,nprint_ene),etot
4559 write (iout,'(//a)') 'Final energies:'
4560 write (iout,'(a4,2x,a12,17a14,3a8)') &
4561 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4562 'RMSnat','NatCONT','NNCONT','RMS'
4566 ist=nres_base(2,ii)+ipatt(2,i)
4568 energia(kk)=ener(kk,ik)
4570 etot=ener(n_ene_comp+1,i)
4571 rmsnat=ener(n_ene_comp+2,i)
4572 rms=ener(n_ene_comp+3,i)
4573 frac=ener(n_ene_comp+4,i)
4574 frac_nn=ener(n_ene_comp+5,i)
4575 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4576 i,str_nam(ii),ist+1,&
4577 (energia(print_order(kk)),kk=1,nprint_ene),&
4578 etot,rmsnat,frac,frac_nn,rms
4580 write (iout,'(/a/)') 'IEXAM array:'
4581 write (iout,'(i5)') nexcl
4583 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4585 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4586 'Max. time for threading step ',max_time_for_thread,&
4587 'Average time for threading step: ',ave_time_for_thread
4589 end subroutine write_thread_summary
4590 !-----------------------------------------------------------------------------
4591 subroutine write_stat_thread(ithread,ipattern,ist)
4593 use energy_data, only: n_ene_comp
4595 ! implicit real*8 (a-h,o-z)
4596 ! include "DIMENSIONS"
4597 ! include "COMMON.CONTROL"
4598 ! include "COMMON.IOUNITS"
4599 ! include "COMMON.THREAD"
4600 ! include "COMMON.FFIELD"
4601 ! include "COMMON.DBASE"
4602 ! include "COMMON.NAMES"
4603 real(kind=8),dimension(0:n_ene) :: energia
4605 integer :: ithread,ipattern,ist,i
4606 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4608 #if defined(AIX) || defined(PGI)
4609 open(istat,file=statname,position='append')
4611 open(istat,file=statname,access='append')
4614 energia(i)=ener(i,ithread)
4616 etot=ener(n_ene_comp+1,ithread)
4617 rmsnat=ener(n_ene_comp+2,ithread)
4618 rms=ener(n_ene_comp+3,ithread)
4619 frac=ener(n_ene_comp+4,ithread)
4620 frac_nn=ener(n_ene_comp+5,ithread)
4621 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4622 ithread,str_nam(ipattern),ist+1,&
4623 (energia(print_order(i)),i=1,nprint_ene),&
4624 etot,rmsnat,frac,frac_nn,rms
4627 end subroutine write_stat_thread
4628 !-----------------------------------------------------------------------------
4630 !-----------------------------------------------------------------------------
4631 end module io_config