8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 integer :: ichir1,ichir2
748 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
749 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
750 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 ! For printing parameters after they are read set the following in the UNRES
755 ! setenv PRINT_PARM YES
757 ! To print parameters in LaTeX format rather than as ASCII tables:
761 call getenv_loc("PRINT_PARM",lancuch)
762 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
763 call getenv_loc("LATEX",lancuch)
764 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
766 dwa16=2.0d0**(1.0d0/6.0d0)
768 ! Assign virtual-bond length
771 vblinv2=vblinv*vblinv
773 ! Read the virtual-bond parameters, masses, and moments of inertia
774 ! and Stokes' radii of the peptide group and side chains
776 allocate(dsc(ntyp1)) !(ntyp1)
777 allocate(dsc_inv(ntyp1)) !(ntyp1)
778 allocate(nbondterm(ntyp)) !(ntyp)
779 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
780 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(msc(ntyp+1)) !(ntyp+1)
782 allocate(isc(ntyp+1)) !(ntyp+1)
783 allocate(restok(ntyp+1)) !(ntyp+1)
784 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(long_r_sidechain(ntyp))
786 allocate(short_r_sidechain(ntyp))
791 read (ibond,*) vbldp0,akp,mp,ip,pstok
794 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
795 dsc(i) = vbldsc0(1,i)
799 dsc_inv(i)=1.0D0/dsc(i)
803 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
805 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
806 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
807 dsc(i) = vbldsc0(1,i)
811 dsc_inv(i)=1.0D0/dsc(i)
816 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
817 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
819 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
821 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
822 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
824 write (iout,'(13x,3f10.5)') &
825 vbldsc0(j,i),aksc(j,i),abond0(j,i)
829 !----------------------------------------------------
830 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
831 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
832 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
833 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
834 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
835 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
845 allocate(liptranene(ntyp))
846 !C reading lipid parameters
847 write (iout,*) "iliptranpar",iliptranpar
849 read(iliptranpar,*) pepliptran
852 read(iliptranpar,*) liptranene(i)
853 print *,liptranene(i)
859 ! Read the parameters of the probability distribution/energy expression
860 ! of the virtual-bond valence angles theta
863 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
864 (bthet(j,i,1,1),j=1,2)
865 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
866 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
867 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
871 athet(1,i,1,-1)=athet(1,i,1,1)
872 athet(2,i,1,-1)=athet(2,i,1,1)
873 bthet(1,i,1,-1)=-bthet(1,i,1,1)
874 bthet(2,i,1,-1)=-bthet(2,i,1,1)
875 athet(1,i,-1,1)=-athet(1,i,1,1)
876 athet(2,i,-1,1)=-athet(2,i,1,1)
877 bthet(1,i,-1,1)=bthet(1,i,1,1)
878 bthet(2,i,-1,1)=bthet(2,i,1,1)
882 athet(1,i,-1,-1)=athet(1,-i,1,1)
883 athet(2,i,-1,-1)=-athet(2,-i,1,1)
884 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
885 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
886 athet(1,i,-1,1)=athet(1,-i,1,1)
887 athet(2,i,-1,1)=-athet(2,-i,1,1)
888 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
889 bthet(2,i,-1,1)=bthet(2,-i,1,1)
890 athet(1,i,1,-1)=-athet(1,-i,1,1)
891 athet(2,i,1,-1)=athet(2,-i,1,1)
892 bthet(1,i,1,-1)=bthet(1,-i,1,1)
893 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
898 polthet(j,i)=polthet(j,-i)
901 gthet(j,i)=gthet(j,-i)
909 'Parameters of the virtual-bond valence angles:'
910 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
911 ' ATHETA0 ',' A1 ',' A2 ',&
914 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
915 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
917 write (iout,'(/a/9x,5a/79(1h-))') &
918 'Parameters of the expression for sigma(theta_c):',&
919 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
920 ' ALPH3 ',' SIGMA0C '
922 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
923 (polthet(j,i),j=0,3),sigc0(i)
925 write (iout,'(/a/9x,5a/79(1h-))') &
926 'Parameters of the second gaussian:',&
927 ' THETA0 ',' SIGMA0 ',' G1 ',&
930 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
931 sig0(i),(gthet(j,i),j=1,3)
935 'Parameters of the virtual-bond valence angles:'
936 write (iout,'(/a/9x,5a/79(1h-))') &
937 'Coefficients of expansion',&
938 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
939 ' b1*10^1 ',' b2*10^1 '
941 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
942 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
943 (10*bthet(j,i,1,1),j=1,2)
945 write (iout,'(/a/9x,5a/79(1h-))') &
946 'Parameters of the expression for sigma(theta_c):',&
947 ' alpha0 ',' alph1 ',' alph2 ',&
948 ' alhp3 ',' sigma0c '
950 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
951 (polthet(j,i),j=0,3),sigc0(i)
953 write (iout,'(/a/9x,5a/79(1h-))') &
954 'Parameters of the second gaussian:',&
955 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
958 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
959 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
965 ! Read the parameters of Utheta determined from ab initio surfaces
966 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
968 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
969 ntheterm3,nsingle,ndouble
970 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
972 !----------------------------------------------------
973 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
974 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
975 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
976 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
977 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
978 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
979 !(maxtheterm,-maxthetyp1:maxthetyp1,&
980 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
981 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
982 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
983 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
984 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
985 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
986 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
987 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
988 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
989 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
990 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
991 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
992 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
993 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
994 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
995 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
996 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
998 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1000 ithetyp(i)=-ithetyp(-i)
1003 aa0thet(:,:,:,:)=0.0d0
1004 aathet(:,:,:,:,:)=0.0d0
1005 bbthet(:,:,:,:,:,:)=0.0d0
1006 ccthet(:,:,:,:,:,:)=0.0d0
1007 ddthet(:,:,:,:,:,:)=0.0d0
1008 eethet(:,:,:,:,:,:)=0.0d0
1009 ffthet(:,:,:,:,:,:,:)=0.0d0
1010 ggthet(:,:,:,:,:,:,:)=0.0d0
1012 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1014 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1015 ! VAR:1=non-glicyne non-proline 2=proline
1016 ! VAR:negative values for D-aminoacid
1018 do j=-nthetyp,nthetyp
1019 do k=-nthetyp,nthetyp
1020 read (ithep,'(6a)',end=111,err=111) res1
1021 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1022 ! VAR: aa0thet is variable describing the average value of Foureir
1023 ! VAR: expansion series
1024 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1025 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1026 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1027 read (ithep,*,end=111,err=111) &
1028 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1029 read (ithep,*,end=111,err=111) &
1030 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1031 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1032 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1033 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1035 read (ithep,*,end=111,err=111) &
1036 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1037 ffthet(lll,llll,ll,i,j,k,iblock),&
1038 ggthet(llll,lll,ll,i,j,k,iblock),&
1039 ggthet(lll,llll,ll,i,j,k,iblock),&
1040 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1045 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1046 ! coefficients of theta-and-gamma-dependent terms are zero.
1047 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1048 ! RECOMENTDED AFTER VERSION 3.3)
1052 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1053 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1055 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1056 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1059 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1061 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1064 ! AND COMMENT THE LOOPS BELOW
1068 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1069 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1071 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1072 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1075 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1077 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1082 ! Substitution for D aminoacids from symmetry.
1085 do j=-nthetyp,nthetyp
1086 do k=-nthetyp,nthetyp
1087 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1089 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1093 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1094 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1095 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1096 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1102 ffthet(llll,lll,ll,i,j,k,iblock)= &
1103 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1104 ffthet(lll,llll,ll,i,j,k,iblock)= &
1105 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1106 ggthet(llll,lll,ll,i,j,k,iblock)= &
1107 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1108 ggthet(lll,llll,ll,i,j,k,iblock)= &
1109 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1118 ! Control printout of the coefficients of virtual-bond-angle potentials
1121 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1126 write (iout,'(//4a)') &
1127 'Type ',onelett(i),onelett(j),onelett(k)
1128 write (iout,'(//a,10x,a)') " l","a[l]"
1129 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1130 write (iout,'(i2,1pe15.5)') &
1131 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1133 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1134 "b",l,"c",l,"d",l,"e",l
1136 write (iout,'(i2,4(1pe15.5))') m,&
1137 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1138 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1142 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1143 "f+",l,"f-",l,"g+",l,"g-",l
1146 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1147 ffthet(n,m,l,i,j,k,iblock),&
1148 ffthet(m,n,l,i,j,k,iblock),&
1149 ggthet(n,m,l,i,j,k,iblock),&
1150 ggthet(m,n,l,i,j,k,iblock)
1160 write (2,*) "Start reading THETA_PDB",ithep_pdb
1162 ! write (2,*) 'i=',i
1163 read (ithep_pdb,*,err=111,end=111) &
1164 a0thet(i),(athet(j,i,1,1),j=1,2),&
1165 (bthet(j,i,1,1),j=1,2)
1166 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1167 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1168 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1169 sigc0(i)=sigc0(i)**2
1172 athet(1,i,1,-1)=athet(1,i,1,1)
1173 athet(2,i,1,-1)=athet(2,i,1,1)
1174 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1175 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1176 athet(1,i,-1,1)=-athet(1,i,1,1)
1177 athet(2,i,-1,1)=-athet(2,i,1,1)
1178 bthet(1,i,-1,1)=bthet(1,i,1,1)
1179 bthet(2,i,-1,1)=bthet(2,i,1,1)
1182 a0thet(i)=a0thet(-i)
1183 athet(1,i,-1,-1)=athet(1,-i,1,1)
1184 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1185 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1186 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1187 athet(1,i,-1,1)=athet(1,-i,1,1)
1188 athet(2,i,-1,1)=-athet(2,-i,1,1)
1189 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1190 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1191 athet(1,i,1,-1)=-athet(1,-i,1,1)
1192 athet(2,i,1,-1)=athet(2,-i,1,1)
1193 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1194 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1195 theta0(i)=theta0(-i)
1199 polthet(j,i)=polthet(j,-i)
1202 gthet(j,i)=gthet(j,-i)
1205 write (2,*) "End reading THETA_PDB"
1210 !-------------------------------------------
1211 allocate(nlob(ntyp1)) !(ntyp1)
1212 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1213 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1214 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1221 gaussc(:,:,:,:)=0.0D0
1225 ! Read the parameters of the probability distribution/energy expression
1226 ! of the side chains.
1229 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1233 dsc_inv(i)=1.0D0/dsc(i)
1244 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1245 ((blower(k,l,1),l=1,k),k=1,3)
1246 censc(1,1,-i)=censc(1,1,i)
1247 censc(2,1,-i)=censc(2,1,i)
1248 censc(3,1,-i)=-censc(3,1,i)
1250 read (irotam,*,end=112,err=112) bsc(j,i)
1251 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1252 ((blower(k,l,j),l=1,k),k=1,3)
1253 censc(1,j,-i)=censc(1,j,i)
1254 censc(2,j,-i)=censc(2,j,i)
1255 censc(3,j,-i)=-censc(3,j,i)
1256 ! BSC is amplitude of Gaussian
1263 akl=akl+blower(k,m,j)*blower(l,m,j)
1267 if (((k.eq.3).and.(l.ne.3)) &
1268 .or.((l.eq.3).and.(k.ne.3))) then
1269 gaussc(k,l,j,-i)=-akl
1270 gaussc(l,k,j,-i)=-akl
1272 gaussc(k,l,j,-i)=akl
1273 gaussc(l,k,j,-i)=akl
1282 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1285 if (nlobi.gt.0) then
1287 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1288 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1289 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1290 'log h',(bsc(j,i),j=1,nlobi)
1291 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1292 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1294 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1295 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1298 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1299 write (iout,'(a,f10.4,4(16x,f10.4))') &
1300 'Center ',(bsc(j,i),j=1,nlobi)
1301 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1310 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1311 ! added by Urszula Kozlowska 07/11/2007
1313 !el Maximum number of SC local term fitting function coefficiants
1314 !el integer,parameter :: maxsccoef=65
1316 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1319 read (irotam,*,end=112,err=112)
1321 read (irotam,*,end=112,err=112)
1324 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1329 ! Read the parameters of the probability distribution/energy expression
1330 ! of the side chains.
1332 write (2,*) "Start reading ROTAM_PDB"
1334 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1338 dsc_inv(i)=1.0D0/dsc(i)
1349 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1350 ((blower(k,l,1),l=1,k),k=1,3)
1352 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1353 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1354 ((blower(k,l,j),l=1,k),k=1,3)
1361 akl=akl+blower(k,m,j)*blower(l,m,j)
1371 write (2,*) "End reading ROTAM_PDB"
1377 ! Read torsional parameters in old format
1379 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1381 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1382 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1383 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1385 !el from energy module--------
1386 allocate(v1(nterm_old,ntortyp,ntortyp))
1387 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1388 !el---------------------------
1393 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1399 write (iout,'(/a/)') 'Torsional constants:'
1402 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1403 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1409 ! Read torsional parameters
1411 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1412 read (itorp,*,end=113,err=113) ntortyp
1413 !el from energy module---------
1414 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1415 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1417 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1418 allocate(vlor2(maxlor,ntortyp,ntortyp))
1419 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1420 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1422 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1423 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1424 !el---------------------------
1427 !el---------------------------
1429 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1431 itortyp(i)=-itortyp(-i)
1433 itortyp(ntyp1)=ntortyp
1434 itortyp(-ntyp1)=-ntortyp
1436 write (iout,*) 'ntortyp',ntortyp
1438 do j=-ntortyp+1,ntortyp-1
1439 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1441 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1442 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1445 do k=1,nterm(i,j,iblock)
1446 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1448 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1449 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1450 v0ij=v0ij+si*v1(k,i,j,iblock)
1452 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1453 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1454 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1456 do k=1,nlor(i,j,iblock)
1457 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1458 vlor2(k,i,j),vlor3(k,i,j)
1459 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1462 v0(-i,-j,iblock)=v0ij
1468 write (iout,'(/a/)') 'Torsional constants:'
1470 do i=-ntortyp,ntortyp
1471 do j=-ntortyp,ntortyp
1472 write (iout,*) 'ityp',i,' jtyp',j
1473 write (iout,*) 'Fourier constants'
1474 do k=1,nterm(i,j,iblock)
1475 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1478 write (iout,*) 'Lorenz constants'
1479 do k=1,nlor(i,j,iblock)
1480 write (iout,'(3(1pe15.5))') &
1481 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1487 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1489 ! 6/23/01 Read parameters for double torsionals
1491 !el from energy module------------
1492 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1493 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1494 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1495 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1496 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1497 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1498 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1499 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1500 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1501 !---------------------------------
1505 do j=-ntortyp+1,ntortyp-1
1506 do k=-ntortyp+1,ntortyp-1
1507 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1508 ! write (iout,*) "OK onelett",
1511 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1512 .or. t3.ne.toronelet(k)) then
1513 write (iout,*) "Error in double torsional parameter file",&
1516 call MPI_Finalize(Ierror)
1518 stop "Error in double torsional parameter file"
1520 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1521 ntermd_2(i,j,k,iblock)
1522 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1523 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1524 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1525 ntermd_1(i,j,k,iblock))
1526 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1527 ntermd_1(i,j,k,iblock))
1528 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1529 ntermd_1(i,j,k,iblock))
1530 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1531 ntermd_1(i,j,k,iblock))
1532 ! Martix of D parameters for one dimesional foureir series
1533 do l=1,ntermd_1(i,j,k,iblock)
1534 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1535 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1536 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1537 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1538 ! write(iout,*) "whcodze" ,
1539 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1541 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1542 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1543 v2s(m,l,i,j,k,iblock),&
1544 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1545 ! Martix of D parameters for two dimesional fourier series
1546 do l=1,ntermd_2(i,j,k,iblock)
1548 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1549 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1550 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1551 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1560 write (iout,*) 'Constants for double torsionals'
1563 do j=-ntortyp+1,ntortyp-1
1564 do k=-ntortyp+1,ntortyp-1
1565 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1566 ' nsingle',ntermd_1(i,j,k,iblock),&
1567 ' ndouble',ntermd_2(i,j,k,iblock)
1569 write (iout,*) 'Single angles:'
1570 do l=1,ntermd_1(i,j,k,iblock)
1571 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1572 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1573 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1574 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1577 write (iout,*) 'Pairs of angles:'
1578 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1579 do l=1,ntermd_2(i,j,k,iblock)
1580 write (iout,'(i5,20f10.5)') &
1581 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1584 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1585 do l=1,ntermd_2(i,j,k,iblock)
1586 write (iout,'(i5,20f10.5)') &
1587 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1588 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1597 ! Read of Side-chain backbone correlation parameters
1598 ! Modified 11 May 2012 by Adasko
1601 read (isccor,*,end=119,err=119) nsccortyp
1603 !el from module energy-------------
1604 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1605 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1606 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1607 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1608 !-----------------------------------
1610 !el from module energy-------------
1611 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1613 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1615 isccortyp(i)=-isccortyp(-i)
1617 iscprol=isccortyp(20)
1618 ! write (iout,*) 'ntortyp',ntortyp
1620 !c maxinter is maximum interaction sites
1621 !el from module energy---------
1622 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1623 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1624 -nsccortyp:nsccortyp))
1625 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1626 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1627 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1628 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1629 !-----------------------------------
1633 read (isccor,*,end=119,err=119) &
1634 nterm_sccor(i,j),nlor_sccor(i,j)
1640 nterm_sccor(-i,j)=nterm_sccor(i,j)
1641 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1642 nterm_sccor(i,-j)=nterm_sccor(i,j)
1643 do k=1,nterm_sccor(i,j)
1644 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1646 if (j.eq.iscprol) then
1647 if (i.eq.isccortyp(10)) then
1648 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1649 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1651 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1652 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1653 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1654 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1655 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1656 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1657 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1658 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1661 if (i.eq.isccortyp(10)) then
1662 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1663 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1665 if (j.eq.isccortyp(10)) then
1666 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1667 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1669 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1670 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1671 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1672 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1673 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1674 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1678 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1679 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1680 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1681 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1684 do k=1,nlor_sccor(i,j)
1685 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1686 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1687 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1688 (1+vlor3sccor(k,i,j)**2)
1690 v0sccor(l,i,j)=v0ijsccor
1691 v0sccor(l,-i,j)=v0ijsccor1
1692 v0sccor(l,i,-j)=v0ijsccor2
1693 v0sccor(l,-i,-j)=v0ijsccor3
1699 !el from module energy-------------
1700 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1702 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1703 ! write (iout,*) 'ntortyp',ntortyp
1705 !c maxinter is maximum interaction sites
1706 !el from module energy---------
1707 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1708 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1709 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1710 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1711 !-----------------------------------
1715 read (isccor,*,end=119,err=119) &
1716 nterm_sccor(i,j),nlor_sccor(i,j)
1720 do k=1,nterm_sccor(i,j)
1721 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1723 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1726 do k=1,nlor_sccor(i,j)
1727 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1728 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1729 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1730 (1+vlor3sccor(k,i,j)**2)
1732 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1740 write (iout,'(/a/)') 'Torsional constants:'
1743 write (iout,*) 'ityp',i,' jtyp',j
1744 write (iout,*) 'Fourier constants'
1745 do k=1,nterm_sccor(i,j)
1746 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1748 write (iout,*) 'Lorenz constants'
1749 do k=1,nlor_sccor(i,j)
1750 write (iout,'(3(1pe15.5))') &
1751 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1758 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1759 ! interaction energy of the Gly, Ala, and Pro prototypes.
1764 write (iout,*) "Coefficients of the cumulants"
1766 read (ifourier,*) nloctyp
1767 !write(iout,*) "nloctyp",nloctyp
1768 !el from module energy-------
1769 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1770 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1771 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1772 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1773 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1774 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1775 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1776 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1780 !--------------------------------
1783 read (ifourier,*,end=115,err=115)
1784 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1786 write (iout,*) 'Type',i
1787 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1797 B1tilde(1,-i) =-b(3)
1799 ! b1tilde(1,i)=0.0d0
1800 ! b1tilde(2,i)=0.0d0
1825 Ctilde(1,2,-i)=-b(9)
1829 ! Ctilde(1,1,i)=0.0d0
1830 ! Ctilde(1,2,i)=0.0d0
1831 ! Ctilde(2,1,i)=0.0d0
1832 ! Ctilde(2,2,i)=0.0d0
1850 Dtilde(1,2,-i)=-b(8)
1854 ! Dtilde(1,1,i)=0.0d0
1855 ! Dtilde(1,2,i)=0.0d0
1856 ! Dtilde(2,1,i)=0.0d0
1857 ! Dtilde(2,2,i)=0.0d0
1858 EE(1,1,i)= b(10)+b(11)
1859 EE(2,2,i)=-b(10)+b(11)
1860 EE(2,1,i)= b(12)-b(13)
1861 EE(1,2,i)= b(12)+b(13)
1862 EE(1,1,-i)= b(10)+b(11)
1863 EE(2,2,-i)=-b(10)+b(11)
1864 EE(2,1,-i)=-b(12)+b(13)
1865 EE(1,2,-i)=-b(12)-b(13)
1871 ! ee(2,1,i)=ee(1,2,i)
1875 write (iout,*) 'Type',i
1877 write(iout,*) B1(1,i),B1(2,i)
1879 write(iout,*) B2(1,i),B2(2,i)
1882 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1886 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1890 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1895 ! Read electrostatic-interaction parameters
1900 write (iout,'(/a)') 'Electrostatic interaction constants:'
1901 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1902 'IT','JT','APP','BPP','AEL6','AEL3'
1904 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1905 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1906 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1907 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1912 app (i,j)=epp(i,j)*rri*rri
1913 bpp (i,j)=-2.0D0*epp(i,j)*rri
1914 ael6(i,j)=elpp6(i,j)*4.2D0**6
1915 ael3(i,j)=elpp3(i,j)*4.2D0**3
1917 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1923 ! Read side-chain interaction parameters.
1925 !el from module energy - COMMON.INTERACT-------
1926 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1927 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1928 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1929 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1930 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1931 allocate(epslip(ntyp,ntyp))
1939 !--------------------------------
1941 read (isidep,*,end=117,err=117) ipot,expon
1942 if (ipot.lt.1 .or. ipot.gt.5) then
1943 write (iout,'(2a)') 'Error while reading SC interaction',&
1944 'potential file - unknown potential type.'
1946 call MPI_Finalize(Ierror)
1952 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1953 ', exponents are ',expon,2*expon
1954 ! goto (10,20,30,30,40) ipot
1956 !----------------------- LJ potential ---------------------------------
1958 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1959 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1960 (sigma0(i),i=1,ntyp)
1962 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1963 write (iout,'(a/)') 'The epsilon array:'
1964 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1965 write (iout,'(/a)') 'One-body parameters:'
1966 write (iout,'(a,4x,a)') 'residue','sigma'
1967 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1970 !----------------------- LJK potential --------------------------------
1972 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1973 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1974 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1976 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1977 write (iout,'(a/)') 'The epsilon array:'
1978 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1979 write (iout,'(/a)') 'One-body parameters:'
1980 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1981 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1985 !---------------------- GB or BP potential -----------------------------
1989 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1991 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1992 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1993 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1994 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1996 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1999 ! For the GB potential convert sigma'**2 into chi'
2002 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2006 write (iout,'(/a/)') 'Parameters of the BP potential:'
2007 write (iout,'(a/)') 'The epsilon array:'
2008 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2009 write (iout,'(/a)') 'One-body parameters:'
2010 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2012 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
2013 chip(i),alp(i),i=1,ntyp)
2016 !--------------------- GBV potential -----------------------------------
2018 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2019 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2020 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2021 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2023 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2024 write (iout,'(a/)') 'The epsilon array:'
2025 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2026 write (iout,'(/a)') 'One-body parameters:'
2027 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2028 's||/s_|_^2',' chip ',' alph '
2029 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2030 sigii(i),chip(i),alp(i),i=1,ntyp)
2033 write(iout,*)"Wrong ipot"
2038 !-----------------------------------------------------------------------
2039 ! Calculate the "working" parameters of SC interactions.
2041 !el from module energy - COMMON.INTERACT-------
2042 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2043 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2044 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2053 !--------------------------------
2058 epslip(i,j)=epslip(j,i)
2063 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2064 sigma(j,i)=sigma(i,j)
2065 rs0(i,j)=dwa16*sigma(i,j)
2069 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2070 'Working parameters of the SC interactions:',&
2071 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2076 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2085 sigeps=dsign(1.0D0,epsij)
2087 aa_aq(i,j)=epsij*rrij*rrij
2088 bb_aq(i,j)=-sigeps*epsij*rrij
2089 aa_aq(j,i)=aa_aq(i,j)
2090 bb_aq(j,i)=bb_aq(i,j)
2091 epsijlip=epslip(i,j)
2092 sigeps=dsign(1.0D0,epsijlip)
2093 epsijlip=dabs(epsijlip)
2094 aa_lip(i,j)=epsijlip*rrij*rrij
2095 bb_lip(i,j)=-sigeps*epsijlip*rrij
2096 aa_lip(j,i)=aa_lip(i,j)
2097 bb_lip(j,i)=bb_lip(i,j)
2098 !C write(iout,*) aa_lip
2100 sigt1sq=sigma0(i)**2
2101 sigt2sq=sigma0(j)**2
2104 ratsig1=sigt2sq/sigt1sq
2105 ratsig2=1.0D0/ratsig1
2106 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2107 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2108 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2112 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2113 sigmaii(i,j)=rsum_max
2114 sigmaii(j,i)=rsum_max
2116 ! sigmaii(i,j)=r0(i,j)
2117 ! sigmaii(j,i)=r0(i,j)
2119 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2120 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2121 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2122 augm(i,j)=epsij*r_augm**(2*expon)
2123 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2130 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2131 restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2132 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2137 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2142 ! Define the SC-p interaction constants (hard-coded; old style)
2145 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2147 ! aad(i,1)=0.3D0*4.0D0**12
2148 ! Following line for constants currently implemented
2149 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2150 aad(i,1)=1.5D0*4.0D0**12
2151 ! aad(i,1)=0.17D0*5.6D0**12
2153 ! "Soft" SC-p repulsion
2155 ! Following line for constants currently implemented
2156 ! aad(i,1)=0.3D0*4.0D0**6
2157 ! "Hard" SC-p repulsion
2158 bad(i,1)=3.0D0*4.0D0**6
2159 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2168 ! 8/9/01 Read the SC-p interaction constants from file
2171 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2174 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2175 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2176 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2177 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2181 write (iout,*) "Parameters of SC-p interactions:"
2183 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2184 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2190 ! Define the constants of the disulfide bridge
2194 ! Old arbitrary potential - commented out.
2199 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2200 ! energy surface of diethyl disulfide.
2201 ! A. Liwo and U. Kozlowska, 11/24/03
2218 write (iout,'(/a)') "Disulfide bridge parameters:"
2219 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2220 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2221 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2222 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2226 111 write (iout,*) "Error reading bending energy parameters."
2228 112 write (iout,*) "Error reading rotamer energy parameters."
2230 113 write (iout,*) "Error reading torsional energy parameters."
2232 114 write (iout,*) "Error reading double torsional energy parameters."
2234 115 write (iout,*) &
2235 "Error reading cumulant (multibody energy) parameters."
2237 116 write (iout,*) "Error reading electrostatic energy parameters."
2239 117 write (iout,*) "Error reading side chain interaction parameters."
2241 118 write (iout,*) "Error reading SCp interaction parameters."
2243 119 write (iout,*) "Error reading SCCOR parameters"
2246 call MPI_Finalize(Ierror)
2250 end subroutine parmread
2252 !-----------------------------------------------------------------------------
2254 !-----------------------------------------------------------------------------
2255 subroutine printmat(ldim,m,n,iout,key,a)
2258 character(len=3),dimension(n) :: key
2259 real(kind=8),dimension(ldim,n) :: a
2261 integer :: i,j,k,m,iout,nlim
2265 write (iout,1000) (key(k),k=i,nlim)
2267 1000 format (/5x,8(6x,a3))
2268 1020 format (/80(1h-)/)
2270 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2273 1010 format (a3,2x,8(f9.4))
2275 end subroutine printmat
2276 !-----------------------------------------------------------------------------
2278 !-----------------------------------------------------------------------------
2280 ! Read the PDB file and convert the peptide geometry into virtual-chain
2283 use energy_data, only: itype
2287 use control, only: rescode
2288 ! implicit real*8 (a-h,o-z)
2289 ! include 'DIMENSIONS'
2290 ! include 'COMMON.LOCAL'
2291 ! include 'COMMON.VAR'
2292 ! include 'COMMON.CHAIN'
2293 ! include 'COMMON.INTERACT'
2294 ! include 'COMMON.IOUNITS'
2295 ! include 'COMMON.GEO'
2296 ! include 'COMMON.NAMES'
2297 ! include 'COMMON.CONTROL'
2298 ! include 'COMMON.DISTFIT'
2299 ! include 'COMMON.SETUP'
2300 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2302 logical :: lprn=.true.,fail
2303 real(kind=8),dimension(3) :: e1,e2,e3
2304 real(kind=8) :: dcj,efree_temp
2305 character(len=3) :: seq,res
2306 character(len=5) :: atom
2307 character(len=80) :: card
2308 real(kind=8),dimension(3,20) :: sccor
2309 integer :: kkk,lll,icha,kupa !rescode,
2312 integer,dimension(2,maxres/3) :: hfrag_alloc
2313 integer,dimension(4,maxres/3) :: bfrag_alloc
2314 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2320 ! write (2,*) "UNRES_PDB",unres_pdb
2328 !-----------------------------
2329 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2330 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2333 read (ipdbin,'(a80)',end=10) card
2334 ! write (iout,'(a)') card
2335 if (card(:5).eq.'HELIX') then
2338 read(card(22:25),*) hfrag(1,nhfrag)
2339 read(card(34:37),*) hfrag(2,nhfrag)
2341 if (card(:5).eq.'SHEET') then
2344 read(card(24:26),*) bfrag(1,nbfrag)
2345 read(card(35:37),*) bfrag(2,nbfrag)
2346 !rc----------------------------------------
2347 !rc to be corrected !!!
2348 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2349 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2350 !rc----------------------------------------
2352 if (card(:3).eq.'END') then
2354 else if (card(:3).eq.'TER') then
2358 itype(ires_old)=ntyp1
2359 itype(ires_old-1)=ntyp1
2361 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2364 dc(j,ires)=sccor(j,iii)
2367 call sccenter(ires,iii,sccor)
2373 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2374 ! Fish out the ATOM cards.
2375 if (index(card(1:4),'ATOM').gt.0) then
2376 read (card(12:16),*) atom
2377 ! write (iout,*) "! ",atom," !",ires
2378 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2379 read (card(23:26),*) ires
2380 read (card(18:20),'(a3)') res
2381 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2382 ! & " ires_old",ires_old
2383 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2384 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2385 if (ires-ishift+ishift1.ne.ires_old) then
2386 ! Calculate the CM of the preceding residue.
2387 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2389 ! write (iout,*) "Calculating sidechain center iii",iii
2392 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2395 call sccenter(ires_old,iii,sccor)
2399 ! Start new residue.
2400 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2403 else if (ibeg.eq.1) then
2404 write (iout,*) "BEG ires",ires
2406 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2410 ires=ires-ishift+ishift1
2412 ! write (iout,*) "ishift",ishift," ires",ires,&
2413 ! " ires_old",ires_old
2415 else if (ibeg.eq.2) then
2417 ishift=-ires_old+ires-1 !!!!!
2418 ishift1=ishift1-1 !!!!!
2419 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2420 ires=ires-ishift+ishift1
2424 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2425 ires=ires-ishift+ishift1
2428 if (res.eq.'ACE' .or. res.eq.'NHE') then
2431 itype(ires)=rescode(ires,res,0)
2434 ires=ires-ishift+ishift1
2436 ! write (iout,*) "ires_old",ires_old," ires",ires
2437 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2440 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2441 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2442 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2443 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2444 ! write (iout,*) "backbone ",atom
2446 write (iout,'(2i3,2x,a,3f8.3)') &
2447 ires,itype(ires),res,(c(j,ires),j=1,3)
2451 sccor(j,iii)=c(j,ires)
2453 ! write (*,*) card(23:27),ires,itype(ires)
2454 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2455 atom.ne.'N' .and. atom.ne.'C' .and. &
2456 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2457 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2458 ! write (iout,*) "sidechain ",atom
2460 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2464 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2465 if (ires.eq.0) return
2466 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2470 ! write (iout,*) i,itype(i)
2471 ! if (itype(i).eq.ntyp1) then
2472 ! write (iout,*) "dummy",i,itype(i)
2474 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2475 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2479 if (itype(i).eq.ntyp1) then
2480 if (itype(i+1).eq.ntyp1) then
2481 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2482 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
2483 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
2485 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2486 ! print *,i,'tu dochodze'
2487 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2495 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2499 dcj=(c(j,i-2)-c(j,i-3))/2.0
2500 if (dcj.eq.0) dcj=1.23591524223
2505 else !itype(i+1).eq.ntyp1
2507 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2508 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2515 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2519 dcj=(c(j,i+3)-c(j,i+2))/2.0
2520 if (dcj.eq.0) dcj=1.23591524223
2525 endif !itype(i+1).eq.ntyp1
2526 endif !itype.eq.ntyp1
2529 ! Calculate the CM of the last side chain.
2533 dc(j,ires)=sccor(j,iii)
2536 call sccenter(ires,iii,sccor)
2542 if (itype(nres).ne.10) then
2546 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2547 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2554 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2558 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2559 c(j,nres)=c(j,nres-1)+dcj
2560 c(j,2*nres)=c(j,nres)
2564 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2566 if (nres.ne.nres0) then
2567 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2569 stop "Error nres value in WHAM input"
2572 !---------------------------------
2573 !el reallocate tables
2576 ! hfrag_alloc(j,i)=hfrag(j,i)
2579 ! bfrag_alloc(j,i)=bfrag(j,i)
2585 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2586 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2587 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2588 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2592 ! hfrag(j,i)=hfrag_alloc(j,i)
2597 ! bfrag(j,i)=bfrag_alloc(j,i)
2600 !el end reallocate tables
2601 !---------------------------------
2609 c(j,2*nres)=c(j,nres)
2611 if (itype(1).eq.ntyp1) then
2615 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2616 call refsys(2,3,4,e1,e2,e3,fail)
2623 c(j,1)=c(j,2)-1.9d0*e2(j)
2627 dcj=(c(j,4)-c(j,3))/2.0
2633 ! Copy the coordinates to reference coordinates
2639 ! Calculate internal coordinates.
2641 write (iout,'(/a)') &
2642 "Cartesian coordinates of the reference structure"
2643 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2644 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2646 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2647 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2648 (c(j,ires+nres),j=1,3)
2651 ! znamy już nres wiec mozna alokowac tablice
2652 ! Calculate internal coordinates.
2653 if(me.eq.king.or..not.out1file)then
2654 write (iout,'(a)') &
2655 "Backbone and SC coordinates as read from the PDB"
2657 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2658 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2659 (c(j,nres+ires),j=1,3)
2663 if(.not.allocated(vbld)) then
2664 allocate(vbld(2*nres))
2669 if(.not.allocated(vbld_inv)) then
2670 allocate(vbld_inv(2*nres))
2676 if(.not.allocated(theta)) then
2677 allocate(theta(nres+2))
2681 if(.not.allocated(phi)) allocate(phi(nres+2))
2682 if(.not.allocated(alph)) allocate(alph(nres+2))
2683 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2684 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2685 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2686 if(.not.allocated(costtab)) allocate(costtab(nres))
2687 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2688 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2689 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2690 if(.not.allocated(xxref)) allocate(xxref(nres))
2691 if(.not.allocated(yyref)) allocate(yyref(nres))
2692 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2693 if(.not.allocated(dc_norm)) then
2694 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2695 allocate(dc_norm(3,0:2*nres+2))
2699 call int_from_cart(.true.,.false.)
2700 call sc_loc_geom(.false.)
2702 thetaref(i)=theta(i)
2712 dc(j,i)=c(j,i+1)-c(j,i)
2713 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2718 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2719 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2721 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2725 ! Copy the coordinates to reference coordinates
2726 ! Splits to single chain if occurs
2730 ! cref(j,i,cou)=c(j,i)
2734 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2735 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2736 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2737 !-----------------------------
2743 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2745 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2748 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2753 cref(j,i,cou)=c(j,i)
2754 cref(j,i+nres,cou)=c(j,i+nres)
2756 chain_rep(j,lll,kkk)=c(j,i)
2757 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2761 write (iout,*) chain_length
2762 if (chain_length.eq.0) chain_length=nres
2764 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2765 chain_rep(j,chain_length+nres,symetr) &
2766 =chain_rep(j,chain_length+nres,1)
2769 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2771 ! do kkk=1,chain_length
2772 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2776 ! makes copy of chains
2777 write (iout,*) "symetr", symetr
2782 if (symetr.gt.1) then
2789 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2795 ! write (iout,*) i,icha
2796 do lll=1,chain_length
2798 if (cou.le.nres) then
2800 kupa=mod(lll,chain_length)
2801 iprzes=(kkk-1)*chain_length+lll
2802 if (kupa.eq.0) kupa=chain_length
2803 ! write (iout,*) "kupa", kupa
2804 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2805 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2812 !-koniec robienia kopii
2815 write (iout,*) "nowa struktura", nperm
2817 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2819 cref(3,i,kkk),cref(1,nres+i,kkk),&
2820 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2822 100 format (//' alpha-carbon coordinates ',&
2823 ' centroid coordinates'/ &
2824 ' ', 6X,'X',11X,'Y',11X,'Z', &
2825 10X,'X',11X,'Y',11X,'Z')
2826 110 format (a,'(',i3,')',6f12.5)
2832 bfrag(i,j)=bfrag(i,j)-ishift
2838 hfrag(i,j)=hfrag(i,j)-ishift
2844 end subroutine readpdb
2845 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2846 !-----------------------------------------------------------------------------
2848 !-----------------------------------------------------------------------------
2849 subroutine read_control
2863 use random, only: random_init
2864 ! implicit real*8 (a-h,o-z)
2865 ! include 'DIMENSIONS'
2867 use prng, only:prng_restart
2869 logical :: OKRandom!, prng_restart
2872 ! include 'COMMON.IOUNITS'
2873 ! include 'COMMON.TIME1'
2874 ! include 'COMMON.THREAD'
2875 ! include 'COMMON.SBRIDGE'
2876 ! include 'COMMON.CONTROL'
2877 ! include 'COMMON.MCM'
2878 ! include 'COMMON.MAP'
2879 ! include 'COMMON.HEADER'
2880 ! include 'COMMON.CSA'
2881 ! include 'COMMON.CHAIN'
2882 ! include 'COMMON.MUCA'
2883 ! include 'COMMON.MD'
2884 ! include 'COMMON.FFIELD'
2885 ! include 'COMMON.INTERACT'
2886 ! include 'COMMON.SETUP'
2887 !el integer :: KDIAG,ICORFL,IXDR
2888 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2889 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2890 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2891 ! character(len=80) :: ucase
2892 character(len=640) :: controlcard
2894 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2900 read (INP,'(a)') titel
2901 call card_concat(controlcard,.true.)
2902 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2903 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2904 call reada(controlcard,'SEED',seed,0.0D0)
2905 call random_init(seed)
2906 ! Set up the time limit (caution! The time must be input in minutes!)
2907 read_cart=index(controlcard,'READ_CART').gt.0
2908 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2909 call readi(controlcard,'SYM',symetr,1)
2910 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2911 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2912 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2913 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2914 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2915 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2916 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2917 call reada(controlcard,'DRMS',drms,0.1D0)
2918 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2919 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2920 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2921 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2922 write (iout,'(a,f10.1)')'DRMS = ',drms
2923 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2924 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2926 call readi(controlcard,'NZ_START',nz_start,0)
2927 call readi(controlcard,'NZ_END',nz_end,0)
2928 call readi(controlcard,'IZ_SC',iz_sc,0)
2929 timlim=60.0D0*timlim
2930 safety = 60.0d0*safety
2933 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2934 !C SHIELD keyword sets if the shielding effect of side-chains is used
2935 !C 0 denots no shielding is used all peptide are equally despite the
2936 !C solvent accesible area
2937 !C 1 the newly introduced function
2938 !C 2 reseved for further possible developement
2939 call readi(controlcard,'SHIELD',shield_mode,0)
2940 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2941 write(iout,*) "shield_mode",shield_mode
2942 !C Varibles set size of box
2943 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2944 call reada(controlcard,'BOXY',boxysize,100.0d0)
2945 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2946 ! CUTOFFF ON ELECTROSTATICS
2947 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2948 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2949 write(iout,*) "R_CUT_ELE=",r_cut_ele
2950 ! Lipidic parameters
2951 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
2952 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
2953 if (lipthick.gt.0.0d0) then
2954 bordliptop=(boxzsize+lipthick)/2.0
2955 bordlipbot=bordliptop-lipthick
2956 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
2957 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
2958 buflipbot=bordlipbot+lipbufthick
2959 bufliptop=bordliptop-lipbufthick
2960 if ((lipbufthick*2.0d0).gt.lipthick) &
2961 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
2962 endif !lipthick.gt.0
2963 write(iout,*) "bordliptop=",bordliptop
2964 write(iout,*) "bordlipbot=",bordlipbot
2965 write(iout,*) "bufliptop=",bufliptop
2966 write(iout,*) "buflipbot=",buflipbot
2967 write (iout,*) "SHIELD MODE",shield_mode
2969 !C-------------------------
2970 minim=(index(controlcard,'MINIMIZE').gt.0)
2971 dccart=(index(controlcard,'CART').gt.0)
2972 overlapsc=(index(controlcard,'OVERLAP').gt.0)
2973 overlapsc=.not.overlapsc
2974 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2975 searchsc=.not.searchsc
2976 sideadd=(index(controlcard,'SIDEADD').gt.0)
2977 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2978 outpdb=(index(controlcard,'PDBOUT').gt.0)
2979 outmol2=(index(controlcard,'MOL2OUT').gt.0)
2980 pdbref=(index(controlcard,'PDBREF').gt.0)
2981 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2982 indpdb=index(controlcard,'PDBSTART')
2983 extconf=(index(controlcard,'EXTCONF').gt.0)
2984 call readi(controlcard,'IPRINT',iprint,0)
2985 call readi(controlcard,'MAXGEN',maxgen,10000)
2986 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2987 call readi(controlcard,"KDIAG",kdiag,0)
2988 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2989 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2990 write (iout,*) "RESCALE_MODE",rescale_mode
2991 split_ene=index(controlcard,'SPLIT_ENE').gt.0
2992 if (index(controlcard,'REGULAR').gt.0.0D0) then
2993 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2997 if (index(controlcard,'CHECKGRAD').gt.0) then
2999 if (index(controlcard,'CART').gt.0) then
3001 elseif (index(controlcard,'CARINT').gt.0) then
3006 elseif (index(controlcard,'THREAD').gt.0) then
3008 call readi(controlcard,'THREAD',nthread,0)
3009 if (nthread.gt.0) then
3010 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3013 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3014 stop 'Error termination in Read_Control.'
3016 else if (index(controlcard,'MCMA').gt.0) then
3018 else if (index(controlcard,'MCEE').gt.0) then
3020 else if (index(controlcard,'MULTCONF').gt.0) then
3022 else if (index(controlcard,'MAP').gt.0) then
3024 call readi(controlcard,'MAP',nmap,0)
3025 else if (index(controlcard,'CSA').gt.0) then
3027 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3029 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3032 !fcm else if (index(controlcard,'MCMF').gt.0) then
3034 else if (index(controlcard,'SOFTREG').gt.0) then
3036 else if (index(controlcard,'CHECK_BOND').gt.0) then
3038 else if (index(controlcard,'TEST').gt.0) then
3040 else if (index(controlcard,'MD').gt.0) then
3042 else if (index(controlcard,'RE ').gt.0) then
3046 lmuca=index(controlcard,'MUCA').gt.0
3047 call readi(controlcard,'MUCADYN',mucadyn,0)
3048 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3049 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3051 write (iout,*) 'MUCADYN=',mucadyn
3052 write (iout,*) 'MUCASMOOTH=',muca_smooth
3055 iscode=index(controlcard,'ONE_LETTER')
3056 indphi=index(controlcard,'PHI')
3057 indback=index(controlcard,'BACK')
3058 iranconf=index(controlcard,'RAND_CONF')
3059 i2ndstr=index(controlcard,'USE_SEC_PRED')
3060 gradout=index(controlcard,'GRADOUT').gt.0
3061 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3062 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3063 if (me.eq.king .or. .not.out1file ) &
3064 write (iout,*) "DISTCHAINMAX",distchainmax
3066 if(me.eq.king.or..not.out1file) &
3067 write (iout,'(2a)') diagmeth(kdiag),&
3068 ' routine used to diagonalize matrices.'
3069 if (shield_mode.gt.0) then
3071 !C VSolvSphere the volume of solving sphere
3073 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3074 !C there will be no distinction between proline peptide group and normal peptide
3075 !C group in case of shielding parameters
3076 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3077 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3078 write (iout,*) VSolvSphere,VSolvSphere_div
3079 !C long axis of side chain
3081 long_r_sidechain(i)=vbldsc0(1,i)
3082 short_r_sidechain(i)=sigma0(i)
3083 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3089 end subroutine read_control
3090 !-----------------------------------------------------------------------------
3091 subroutine read_REMDpar
3093 ! Read REMD settings
3100 use control_data, only:out1file
3101 ! implicit real*8 (a-h,o-z)
3102 ! include 'DIMENSIONS'
3103 ! include 'COMMON.IOUNITS'
3104 ! include 'COMMON.TIME1'
3105 ! include 'COMMON.MD'
3108 !el include 'COMMON.LANGEVIN'
3110 !el include 'COMMON.LANGEVIN.lang0'
3112 ! include 'COMMON.INTERACT'
3113 ! include 'COMMON.NAMES'
3114 ! include 'COMMON.GEO'
3115 ! include 'COMMON.REMD'
3116 ! include 'COMMON.CONTROL'
3117 ! include 'COMMON.SETUP'
3118 ! character(len=80) :: ucase
3119 character(len=320) :: controlcard
3120 character(len=3200) :: controlcard1
3121 integer :: iremd_m_total
3124 ! real(kind=8) :: var,ene
3126 if(me.eq.king.or..not.out1file) &
3127 write (iout,*) "REMD setup"
3129 call card_concat(controlcard,.true.)
3130 call readi(controlcard,"NREP",nrep,3)
3131 call readi(controlcard,"NSTEX",nstex,1000)
3132 call reada(controlcard,"RETMIN",retmin,10.0d0)
3133 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3134 mremdsync=(index(controlcard,'SYNC').gt.0)
3135 call readi(controlcard,"NSYN",i_sync_step,100)
3136 restart1file=(index(controlcard,'REST1FILE').gt.0)
3137 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3138 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3139 if(max_cache_traj_use.gt.max_cache_traj) &
3140 max_cache_traj_use=max_cache_traj
3141 if(me.eq.king.or..not.out1file) then
3142 !d if (traj1file) then
3143 !rc caching is in testing - NTWX is not ignored
3144 !d write (iout,*) "NTWX value is ignored"
3145 !d write (iout,*) " trajectory is stored to one file by master"
3146 !d write (iout,*) " before exchange at NSTEX intervals"
3148 write (iout,*) "NREP= ",nrep
3149 write (iout,*) "NSTEX= ",nstex
3150 write (iout,*) "SYNC= ",mremdsync
3151 write (iout,*) "NSYN= ",i_sync_step
3152 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3155 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3156 if (index(controlcard,'TLIST').gt.0) then
3158 call card_concat(controlcard1,.true.)
3159 read(controlcard1,*) (remd_t(i),i=1,nrep)
3160 if(me.eq.king.or..not.out1file) &
3161 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3164 if (index(controlcard,'MLIST').gt.0) then
3166 call card_concat(controlcard1,.true.)
3167 read(controlcard1,*) (remd_m(i),i=1,nrep)
3168 if(me.eq.king.or..not.out1file) then
3169 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3172 iremd_m_total=iremd_m_total+remd_m(i)
3174 write (iout,*) 'Total number of replicas ',iremd_m_total
3177 if(me.eq.king.or..not.out1file) &
3178 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3180 end subroutine read_REMDpar
3181 !-----------------------------------------------------------------------------
3182 subroutine read_MDpar
3186 use control_data, only: r_cut,rlamb,out1file
3188 use geometry_data, only: pi
3190 ! implicit real*8 (a-h,o-z)
3191 ! include 'DIMENSIONS'
3192 ! include 'COMMON.IOUNITS'
3193 ! include 'COMMON.TIME1'
3194 ! include 'COMMON.MD'
3197 !el include 'COMMON.LANGEVIN'
3199 !el include 'COMMON.LANGEVIN.lang0'
3201 ! include 'COMMON.INTERACT'
3202 ! include 'COMMON.NAMES'
3203 ! include 'COMMON.GEO'
3204 ! include 'COMMON.SETUP'
3205 ! include 'COMMON.CONTROL'
3206 ! include 'COMMON.SPLITELE'
3207 ! character(len=80) :: ucase
3208 character(len=320) :: controlcard
3213 call card_concat(controlcard,.true.)
3214 call readi(controlcard,"NSTEP",n_timestep,1000000)
3215 call readi(controlcard,"NTWE",ntwe,100)
3216 call readi(controlcard,"NTWX",ntwx,1000)
3217 call reada(controlcard,"DT",d_time,1.0d-1)
3218 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3219 call reada(controlcard,"DAMAX",damax,1.0d1)
3220 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3221 call readi(controlcard,"LANG",lang,0)
3222 RESPA = index(controlcard,"RESPA") .gt. 0
3223 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3224 ntime_split0=ntime_split
3225 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3226 ntime_split0=ntime_split
3227 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3228 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3229 rest = index(controlcard,"REST").gt.0
3230 tbf = index(controlcard,"TBF").gt.0
3231 usampl = index(controlcard,"USAMPL").gt.0
3232 mdpdb = index(controlcard,"MDPDB").gt.0
3233 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3234 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3235 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3236 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3237 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3238 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3239 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3240 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3241 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3242 large = index(controlcard,"LARGE").gt.0
3243 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3244 rattle = index(controlcard,"RATTLE").gt.0
3245 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3251 if(me.eq.king.or..not.out1file) then
3253 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3255 write (iout,'(a)') "The units are:"
3256 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3257 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3258 " acceleration: angstrom/(48.9 fs)**2"
3259 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3261 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3262 write (iout,'(a60,f10.5,a)') &
3263 "Initial time step of numerical integration:",d_time,&
3265 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3267 write (iout,'(2a,i4,a)') &
3268 "A-MTS algorithm used; initial time step for fast-varying",&
3269 " short-range forces split into",ntime_split," steps."
3270 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3271 r_cut," lambda",rlamb
3273 write (iout,'(2a,f10.5)') &
3274 "Maximum acceleration threshold to reduce the time step",&
3275 "/increase split number:",damax
3276 write (iout,'(2a,f10.5)') &
3277 "Maximum predicted energy drift to reduce the timestep",&
3278 "/increase split number:",edriftmax
3279 write (iout,'(a60,f10.5)') &
3280 "Maximum velocity threshold to reduce velocities:",dvmax
3281 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3282 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3283 if (rattle) write (iout,'(a60)') &
3284 "Rattle algorithm used to constrain the virtual bonds"
3288 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3289 call reada(controlcard,"RWAT",rwat,1.4d0)
3290 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3291 surfarea=index(controlcard,"SURFAREA").gt.0
3292 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3293 if(me.eq.king.or..not.out1file)then
3294 write (iout,'(/a,$)') "Langevin dynamics calculation"
3296 write (iout,'(a/)') &
3297 " with direct integration of Langevin equations"
3298 else if (lang.eq.2) then
3299 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3300 else if (lang.eq.3) then
3301 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3302 else if (lang.eq.4) then
3303 write (iout,'(a/)') " in overdamped mode"
3305 write (iout,'(//a,i5)') &
3306 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3309 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3310 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3311 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3312 write (iout,'(a60,f10.5)') &
3313 "Scaling factor of the friction forces:",scal_fric
3314 if (surfarea) write (iout,'(2a,i10,a)') &
3315 "Friction coefficients will be scaled by solvent-accessible",&
3316 " surface area every",reset_fricmat," steps."
3318 ! Calculate friction coefficients and bounds of stochastic forces
3319 eta=6*pi*cPoise*etawat
3320 if(me.eq.king.or..not.out1file) &
3321 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3323 gamp=scal_fric*(pstok+rwat)*eta
3324 stdfp=dsqrt(2*Rb*t_bath/d_time)
3325 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3327 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3328 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3330 if(me.eq.king.or..not.out1file)then
3331 write (iout,'(/2a/)') &
3332 "Radii of site types and friction coefficients and std's of",&
3333 " stochastic forces of fully exposed sites"
3334 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3336 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3337 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3341 if(me.eq.king.or..not.out1file)then
3342 write (iout,'(a)') "Berendsen bath calculation"
3343 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3344 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3346 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3347 count_reset_moment," steps"
3349 write (iout,'(a,i10,a)') &
3350 "Velocities will be reset at random every",count_reset_vel,&
3354 if(me.eq.king.or..not.out1file) &
3355 write (iout,'(a31)') "Microcanonical mode calculation"
3357 if(me.eq.king.or..not.out1file)then
3358 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3360 write(iout,*) "MD running with constraints."
3361 write(iout,*) "Equilibration time ", eq_time, " mtus."
3362 write(iout,*) "Constraining ", nfrag," fragments."
3363 write(iout,*) "Length of each fragment, weight and q0:"
3365 write (iout,*) "Set of restraints #",iset
3367 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3368 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3370 write(iout,*) "constraints between ", npair, "fragments."
3371 write(iout,*) "constraint pairs, weights and q0:"
3373 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3374 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3376 write(iout,*) "angle constraints within ", nfrag_back,&
3377 "backbone fragments."
3378 write(iout,*) "fragment, weights:"
3380 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3381 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3382 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3385 iset=mod(kolor,nset)+1
3388 if(me.eq.king.or..not.out1file) &
3389 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3391 end subroutine read_MDpar
3392 !-----------------------------------------------------------------------------
3396 ! implicit real*8 (a-h,o-z)
3397 ! include 'DIMENSIONS'
3398 ! include 'COMMON.MAP'
3399 ! include 'COMMON.IOUNITS'
3400 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3401 character(len=80) :: mapcard !,ucase
3404 ! real(kind=8) :: var,ene
3407 read (inp,'(a)') mapcard
3408 mapcard=ucase(mapcard)
3409 if (index(mapcard,'PHI').gt.0) then
3411 else if (index(mapcard,'THE').gt.0) then
3413 else if (index(mapcard,'ALP').gt.0) then
3415 else if (index(mapcard,'OME').gt.0) then
3418 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3419 stop 'Error - illegal variable spec in MAP card.'
3421 call readi (mapcard,'RES1',res1(imap),0)
3422 call readi (mapcard,'RES2',res2(imap),0)
3423 if (res1(imap).eq.0) then
3424 res1(imap)=res2(imap)
3425 else if (res2(imap).eq.0) then
3426 res2(imap)=res1(imap)
3428 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3429 write (iout,'(a)') &
3430 'Error - illegal definition of variable group in MAP.'
3431 stop 'Error - illegal definition of variable group in MAP.'
3433 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3434 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3435 call readi(mapcard,'NSTEP',nstep(imap),0)
3436 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3437 write (iout,'(a)') &
3438 'Illegal boundary and/or step size specification in MAP.'
3439 stop 'Illegal boundary and/or step size specification in MAP.'
3443 end subroutine map_read
3444 !-----------------------------------------------------------------------------
3447 use control_data, only: vdisulf
3449 ! implicit real*8 (a-h,o-z)
3450 ! include 'DIMENSIONS'
3451 ! include 'COMMON.IOUNITS'
3452 ! include 'COMMON.GEO'
3453 ! include 'COMMON.CSA'
3454 ! include 'COMMON.BANK'
3455 ! include 'COMMON.CONTROL'
3456 ! character(len=80) :: ucase
3457 character(len=620) :: mcmcard
3459 ! integer :: ntf,ik,iw_pdb
3460 ! real(kind=8) :: var,ene
3462 call card_concat(mcmcard,.true.)
3464 call readi(mcmcard,'NCONF',nconf,50)
3465 call readi(mcmcard,'NADD',nadd,0)
3466 call readi(mcmcard,'JSTART',jstart,1)
3467 call readi(mcmcard,'JEND',jend,1)
3468 call readi(mcmcard,'NSTMAX',nstmax,500000)
3469 call readi(mcmcard,'N0',n0,1)
3470 call readi(mcmcard,'N1',n1,6)
3471 call readi(mcmcard,'N2',n2,4)
3472 call readi(mcmcard,'N3',n3,0)
3473 call readi(mcmcard,'N4',n4,0)
3474 call readi(mcmcard,'N5',n5,0)
3475 call readi(mcmcard,'N6',n6,10)
3476 call readi(mcmcard,'N7',n7,0)
3477 call readi(mcmcard,'N8',n8,0)
3478 call readi(mcmcard,'N9',n9,0)
3479 call readi(mcmcard,'N14',n14,0)
3480 call readi(mcmcard,'N15',n15,0)
3481 call readi(mcmcard,'N16',n16,0)
3482 call readi(mcmcard,'N17',n17,0)
3483 call readi(mcmcard,'N18',n18,0)
3485 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3487 call readi(mcmcard,'NDIFF',ndiff,2)
3488 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3489 call readi(mcmcard,'IS1',is1,1)
3490 call readi(mcmcard,'IS2',is2,8)
3491 call readi(mcmcard,'NRAN0',nran0,4)
3492 call readi(mcmcard,'NRAN1',nran1,2)
3493 call readi(mcmcard,'IRR',irr,1)
3494 call readi(mcmcard,'NSEED',nseed,20)
3495 call readi(mcmcard,'NTOTAL',ntotal,10000)
3496 call reada(mcmcard,'CUT1',cut1,2.0d0)
3497 call reada(mcmcard,'CUT2',cut2,5.0d0)
3498 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3499 call readi(mcmcard,'ICMAX',icmax,3)
3500 call readi(mcmcard,'IRESTART',irestart,0)
3501 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3504 call reada(mcmcard,'DELE',dele,20.0d0)
3505 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3506 call readi(mcmcard,'IREF',iref,0)
3507 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3508 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3509 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3510 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3511 write (iout,*) "NCONF_IN",nconf_in
3513 end subroutine csaread
3514 !-----------------------------------------------------------------------------
3518 use control_data, only: MaxMoveType
3521 ! implicit real*8 (a-h,o-z)
3522 ! include 'DIMENSIONS'
3523 ! include 'COMMON.MCM'
3524 ! include 'COMMON.MCE'
3525 ! include 'COMMON.IOUNITS'
3526 ! character(len=80) :: ucase
3527 character(len=320) :: mcmcard
3530 ! real(kind=8) :: var,ene
3532 call card_concat(mcmcard,.true.)
3533 call readi(mcmcard,'MAXACC',maxacc,100)
3534 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3535 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3536 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3537 call readi(mcmcard,'MAXREPM',maxrepm,200)
3538 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3539 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3540 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3541 call reada(mcmcard,'E_UP',e_up,5.0D0)
3542 call reada(mcmcard,'DELTE',delte,0.1D0)
3543 call readi(mcmcard,'NSWEEP',nsweep,5)
3544 call readi(mcmcard,'NSTEPH',nsteph,0)
3545 call readi(mcmcard,'NSTEPC',nstepc,0)
3546 call reada(mcmcard,'TMIN',tmin,298.0D0)
3547 call reada(mcmcard,'TMAX',tmax,298.0D0)
3548 call readi(mcmcard,'NWINDOW',nwindow,0)
3549 call readi(mcmcard,'PRINT_MC',print_mc,0)
3550 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3551 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3552 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3553 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3554 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3555 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3556 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3557 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3558 if (nwindow.gt.0) then
3559 allocate(winstart(nwindow)) !!el (maxres)
3560 allocate(winend(nwindow)) !!el
3561 allocate(winlen(nwindow)) !!el
3562 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3564 winlen(i)=winend(i)-winstart(i)+1
3567 if (tmax.lt.tmin) tmax=tmin
3568 if (tmax.eq.tmin) then
3572 if (nstepc.gt.0 .and. nsteph.gt.0) then
3573 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3574 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3576 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3577 ! Probabilities of different move types
3578 sumpro_type(0)=0.0D0
3579 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3580 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3581 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3582 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3583 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3584 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3585 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3587 print *,'i',i,' sumprotype',sumpro_type(i)
3588 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3589 print *,'i',i,' sumprotype',sumpro_type(i)
3592 end subroutine mcmread
3593 !-----------------------------------------------------------------------------
3594 subroutine read_minim
3597 ! implicit real*8 (a-h,o-z)
3598 ! include 'DIMENSIONS'
3599 ! include 'COMMON.MINIM'
3600 ! include 'COMMON.IOUNITS'
3601 ! character(len=80) :: ucase
3602 character(len=320) :: minimcard
3604 ! integer :: ntf,ik,iw_pdb
3605 ! real(kind=8) :: var,ene
3607 call card_concat(minimcard,.true.)
3608 call readi(minimcard,'MAXMIN',maxmin,2000)
3609 call readi(minimcard,'MAXFUN',maxfun,5000)
3610 call readi(minimcard,'MINMIN',minmin,maxmin)
3611 call readi(minimcard,'MINFUN',minfun,maxmin)
3612 call reada(minimcard,'TOLF',tolf,1.0D-2)
3613 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3614 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3615 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3616 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3617 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3618 'Options in energy minimization:'
3619 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3620 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3621 'MinMin:',MinMin,' MinFun:',MinFun,&
3622 ' TolF:',TolF,' RTolF:',RTolF
3624 end subroutine read_minim
3625 !-----------------------------------------------------------------------------
3626 subroutine openunits
3628 use MD_data, only: usampl
3631 use control_data, only:out1file
3632 use control, only: getenv_loc
3633 ! implicit real*8 (a-h,o-z)
3634 ! include 'DIMENSIONS'
3637 character(len=16) :: form,nodename
3638 integer :: nodelen,ierror,npos
3640 ! include 'COMMON.SETUP'
3641 ! include 'COMMON.IOUNITS'
3642 ! include 'COMMON.MD'
3643 ! include 'COMMON.CONTROL'
3644 integer :: lenpre,lenpot,lentmp !,ilen
3646 character(len=3) :: out1file_text !,ucase
3647 character(len=3) :: ll
3650 ! integer :: ntf,ik,iw_pdb
3651 ! real(kind=8) :: var,ene
3653 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3654 call getenv_loc("PREFIX",prefix)
3656 call getenv_loc("POT",pot)
3657 call getenv_loc("DIRTMP",tmpdir)
3658 call getenv_loc("CURDIR",curdir)
3659 call getenv_loc("OUT1FILE",out1file_text)
3660 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3661 out1file_text=ucase(out1file_text)
3662 if (out1file_text(1:1).eq."Y") then
3665 out1file=fg_rank.gt.0
3670 if (lentmp.gt.0) then
3671 write (*,'(80(1h!))')
3672 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3673 write (*,'(80(1h!))')
3674 write (*,*)"All output files will be on node /tmp directory."
3676 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3677 if (me.eq.king) then
3678 write (*,*) "The master node is ",nodename
3679 else if (fg_rank.eq.0) then
3680 write (*,*) "I am the CG slave node ",nodename
3682 write (*,*) "I am the FG slave node ",nodename
3685 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3686 lenpre = lentmp+lenpre+1
3688 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3689 ! Get the names and open the input files
3690 #if defined(WINIFL) || defined(WINPGI)
3691 open(1,file=pref_orig(:ilen(pref_orig))// &
3692 '.inp',status='old',readonly,shared)
3693 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3694 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3695 ! Get parameter filenames and open the parameter files.
3696 call getenv_loc('BONDPAR',bondname)
3697 open (ibond,file=bondname,status='old',readonly,shared)
3698 call getenv_loc('THETPAR',thetname)
3699 open (ithep,file=thetname,status='old',readonly,shared)
3700 call getenv_loc('ROTPAR',rotname)
3701 open (irotam,file=rotname,status='old',readonly,shared)
3702 call getenv_loc('TORPAR',torname)
3703 open (itorp,file=torname,status='old',readonly,shared)
3704 call getenv_loc('TORDPAR',tordname)
3705 open (itordp,file=tordname,status='old',readonly,shared)
3706 call getenv_loc('FOURIER',fouriername)
3707 open (ifourier,file=fouriername,status='old',readonly,shared)
3708 call getenv_loc('ELEPAR',elename)
3709 open (ielep,file=elename,status='old',readonly,shared)
3710 call getenv_loc('SIDEPAR',sidename)
3711 open (isidep,file=sidename,status='old',readonly,shared)
3712 #elif (defined CRAY) || (defined AIX)
3713 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3715 ! print *,"Processor",myrank," opened file 1"
3716 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3717 ! print *,"Processor",myrank," opened file 9"
3718 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3719 ! Get parameter filenames and open the parameter files.
3720 call getenv_loc('BONDPAR',bondname)
3721 open (ibond,file=bondname,status='old',action='read')
3722 ! print *,"Processor",myrank," opened file IBOND"
3723 call getenv_loc('THETPAR',thetname)
3724 open (ithep,file=thetname,status='old',action='read')
3725 ! print *,"Processor",myrank," opened file ITHEP"
3726 call getenv_loc('ROTPAR',rotname)
3727 open (irotam,file=rotname,status='old',action='read')
3728 ! print *,"Processor",myrank," opened file IROTAM"
3729 call getenv_loc('TORPAR',torname)
3730 open (itorp,file=torname,status='old',action='read')
3731 ! print *,"Processor",myrank," opened file ITORP"
3732 call getenv_loc('TORDPAR',tordname)
3733 open (itordp,file=tordname,status='old',action='read')
3734 ! print *,"Processor",myrank," opened file ITORDP"
3735 call getenv_loc('SCCORPAR',sccorname)
3736 open (isccor,file=sccorname,status='old',action='read')
3737 ! print *,"Processor",myrank," opened file ISCCOR"
3738 call getenv_loc('FOURIER',fouriername)
3739 open (ifourier,file=fouriername,status='old',action='read')
3740 ! print *,"Processor",myrank," opened file IFOURIER"
3741 call getenv_loc('ELEPAR',elename)
3742 open (ielep,file=elename,status='old',action='read')
3743 ! print *,"Processor",myrank," opened file IELEP"
3744 call getenv_loc('SIDEPAR',sidename)
3745 open (isidep,file=sidename,status='old',action='read')
3746 call getenv_loc('LIPTRANPAR',liptranname)
3747 open (iliptranpar,file=liptranname,status='old',action='read')
3748 ! print *,"Processor",myrank," opened file ISIDEP"
3749 ! print *,"Processor",myrank," opened parameter files"
3751 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3752 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3753 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3754 ! Get parameter filenames and open the parameter files.
3755 call getenv_loc('BONDPAR',bondname)
3756 open (ibond,file=bondname,status='old')
3757 call getenv_loc('THETPAR',thetname)
3758 open (ithep,file=thetname,status='old')
3759 call getenv_loc('ROTPAR',rotname)
3760 open (irotam,file=rotname,status='old')
3761 call getenv_loc('TORPAR',torname)
3762 open (itorp,file=torname,status='old')
3763 call getenv_loc('TORDPAR',tordname)
3764 open (itordp,file=tordname,status='old')
3765 call getenv_loc('SCCORPAR',sccorname)
3766 open (isccor,file=sccorname,status='old')
3767 call getenv_loc('FOURIER',fouriername)
3768 open (ifourier,file=fouriername,status='old')
3769 call getenv_loc('ELEPAR',elename)
3770 open (ielep,file=elename,status='old')
3771 call getenv_loc('SIDEPAR',sidename)
3772 open (isidep,file=sidename,status='old')
3773 call getenv_loc('LIPTRANPAR',liptranname)
3774 open (iliptranpar,file=liptranname,status='old')
3776 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3778 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3779 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3780 ! Get parameter filenames and open the parameter files.
3781 call getenv_loc('BONDPAR',bondname)
3782 open (ibond,file=bondname,status='old',action='read')
3783 call getenv_loc('THETPAR',thetname)
3784 open (ithep,file=thetname,status='old',action='read')
3785 call getenv_loc('ROTPAR',rotname)
3786 open (irotam,file=rotname,status='old',action='read')
3787 call getenv_loc('TORPAR',torname)
3788 open (itorp,file=torname,status='old',action='read')
3789 call getenv_loc('TORDPAR',tordname)
3790 open (itordp,file=tordname,status='old',action='read')
3791 call getenv_loc('SCCORPAR',sccorname)
3792 open (isccor,file=sccorname,status='old',action='read')
3794 call getenv_loc('THETPARPDB',thetname_pdb)
3795 print *,"thetname_pdb ",thetname_pdb
3796 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3797 print *,ithep_pdb," opened"
3799 call getenv_loc('FOURIER',fouriername)
3800 open (ifourier,file=fouriername,status='old',readonly)
3801 call getenv_loc('ELEPAR',elename)
3802 open (ielep,file=elename,status='old',readonly)
3803 call getenv_loc('SIDEPAR',sidename)
3804 open (isidep,file=sidename,status='old',readonly)
3805 call getenv_loc('LIPTRANPAR',liptranname)
3806 open (iliptranpar,file=liptranname,status='old',action='read')
3808 call getenv_loc('ROTPARPDB',rotname_pdb)
3809 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3814 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3815 ! Use -DOLDSCP to use hard-coded constants instead.
3817 call getenv_loc('SCPPAR',scpname)
3818 #if defined(WINIFL) || defined(WINPGI)
3819 open (iscpp,file=scpname,status='old',readonly,shared)
3820 #elif (defined CRAY) || (defined AIX)
3821 open (iscpp,file=scpname,status='old',action='read')
3823 open (iscpp,file=scpname,status='old')
3825 open (iscpp,file=scpname,status='old',action='read')
3828 call getenv_loc('PATTERN',patname)
3829 #if defined(WINIFL) || defined(WINPGI)
3830 open (icbase,file=patname,status='old',readonly,shared)
3831 #elif (defined CRAY) || (defined AIX)
3832 open (icbase,file=patname,status='old',action='read')
3834 open (icbase,file=patname,status='old')
3836 open (icbase,file=patname,status='old',action='read')
3839 ! Open output file only for CG processes
3840 ! print *,"Processor",myrank," fg_rank",fg_rank
3841 if (fg_rank.eq.0) then
3843 if (nodes.eq.1) then
3846 npos = dlog10(dfloat(nodes-1))+1
3848 if (npos.lt.3) npos=3
3849 write (liczba,'(i1)') npos
3850 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3852 write (liczba,form) me
3853 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3854 liczba(:ilen(liczba))
3855 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3857 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3859 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3860 liczba(:ilen(liczba))//'.mol2'
3861 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3862 liczba(:ilen(liczba))//'.stat'
3864 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3865 //liczba(:ilen(liczba))//'.stat')
3866 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3869 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3870 liczba(:ilen(liczba))//'.const'
3875 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3876 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3877 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3878 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3879 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3881 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3883 rest2name=prefix(:ilen(prefix))//'.rst'
3885 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3888 #if defined(AIX) || defined(PGI)
3889 if (me.eq.king .or. .not. out1file) &
3890 open(iout,file=outname,status='unknown')
3892 if (fg_rank.gt.0) then
3893 write (liczba,'(i3.3)') myrank/nfgtasks
3894 write (ll,'(bz,i3.3)') fg_rank
3895 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3900 open(igeom,file=intname,status='unknown',position='append')
3901 open(ipdb,file=pdbname,status='unknown')
3902 open(imol2,file=mol2name,status='unknown')
3903 open(istat,file=statname,status='unknown',position='append')
3905 !1out open(iout,file=outname,status='unknown')
3908 if (me.eq.king .or. .not.out1file) &
3909 open(iout,file=outname,status='unknown')
3911 if (fg_rank.gt.0) then
3912 write (liczba,'(i3.3)') myrank/nfgtasks
3913 write (ll,'(bz,i3.3)') fg_rank
3914 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3919 open(igeom,file=intname,status='unknown',access='append')
3920 open(ipdb,file=pdbname,status='unknown')
3921 open(imol2,file=mol2name,status='unknown')
3922 open(istat,file=statname,status='unknown',access='append')
3924 !1out open(iout,file=outname,status='unknown')
3927 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3928 csa_seed=prefix(:lenpre)//'.CSA.seed'
3929 csa_history=prefix(:lenpre)//'.CSA.history'
3930 csa_bank=prefix(:lenpre)//'.CSA.bank'
3931 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3932 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3933 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3934 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3935 csa_int=prefix(:lenpre)//'.int'
3936 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3937 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3938 csa_in=prefix(:lenpre)//'.CSA.in'
3939 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3942 write (iout,'(80(1h-))')
3943 write (iout,'(30x,a)') "FILE ASSIGNMENT"
3944 write (iout,'(80(1h-))')
3945 write (iout,*) "Input file : ",&
3946 pref_orig(:ilen(pref_orig))//'.inp'
3947 write (iout,*) "Output file : ",&
3948 outname(:ilen(outname))
3950 write (iout,*) "Sidechain potential file : ",&
3951 sidename(:ilen(sidename))
3953 write (iout,*) "SCp potential file : ",&
3954 scpname(:ilen(scpname))
3956 write (iout,*) "Electrostatic potential file : ",&
3957 elename(:ilen(elename))
3958 write (iout,*) "Cumulant coefficient file : ",&
3959 fouriername(:ilen(fouriername))
3960 write (iout,*) "Torsional parameter file : ",&
3961 torname(:ilen(torname))
3962 write (iout,*) "Double torsional parameter file : ",&
3963 tordname(:ilen(tordname))
3964 write (iout,*) "SCCOR parameter file : ",&
3965 sccorname(:ilen(sccorname))
3966 write (iout,*) "Bond & inertia constant file : ",&
3967 bondname(:ilen(bondname))
3968 write (iout,*) "Bending parameter file : ",&
3969 thetname(:ilen(thetname))
3970 write (iout,*) "Rotamer parameter file : ",&
3971 rotname(:ilen(rotname))
3974 write (iout,*) "Thetpdb parameter file : ",&
3975 thetname_pdb(:ilen(thetname_pdb))
3978 write (iout,*) "Threading database : ",&
3979 patname(:ilen(patname))
3981 write (iout,*)" DIRTMP : ",&
3983 write (iout,'(80(1h-))')
3986 end subroutine openunits
3987 !-----------------------------------------------------------------------------
3990 use geometry_data, only: nres,dc
3992 ! implicit real*8 (a-h,o-z)
3993 ! include 'DIMENSIONS'
3994 ! include 'COMMON.CHAIN'
3995 ! include 'COMMON.IOUNITS'
3996 ! include 'COMMON.MD'
3999 ! real(kind=8) :: var,ene
4001 open(irest2,file=rest2name,status='unknown')
4002 read(irest2,*) totT,EK,potE,totE,t_bath
4004 ! AL 4/17/17: Now reading d_t(0,:) too
4006 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4009 ! AL 4/17/17: Now reading d_c(0,:) too
4011 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4014 read (irest2,*) iset
4018 end subroutine readrst
4019 !-----------------------------------------------------------------------------
4020 subroutine read_fragments
4024 use control_data, only:out1file
4027 ! implicit real*8 (a-h,o-z)
4028 ! include 'DIMENSIONS'
4032 ! include 'COMMON.SETUP'
4033 ! include 'COMMON.CHAIN'
4034 ! include 'COMMON.IOUNITS'
4035 ! include 'COMMON.MD'
4036 ! include 'COMMON.CONTROL'
4039 ! real(kind=8) :: var,ene
4041 read(inp,*) nset,nfrag,npair,nfrag_back
4043 !el from module energy
4044 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4045 if(.not.allocated(wfrag_back)) then
4046 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4047 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4049 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4050 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4052 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4053 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4056 if(me.eq.king.or..not.out1file) &
4057 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4058 " nfrag_back",nfrag_back
4060 read(inp,*) mset(iset)
4062 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4064 if(me.eq.king.or..not.out1file) &
4065 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4066 ifrag(2,i,iset), qinfrag(i,iset)
4069 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4071 if(me.eq.king.or..not.out1file) &
4072 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4073 ipair(2,i,iset), qinpair(i,iset)
4076 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4077 wfrag_back(3,i,iset),&
4078 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4079 if(me.eq.king.or..not.out1file) &
4080 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4081 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4085 end subroutine read_fragments
4086 !-----------------------------------------------------------------------------
4088 !-----------------------------------------------------------------------------
4092 ! implicit real*8 (a-h,o-z)
4093 ! include 'DIMENSIONS'
4094 ! include 'COMMON.CSA'
4095 ! include 'COMMON.BANK'
4096 ! include 'COMMON.IOUNITS'
4098 ! integer :: ntf,ik,iw_pdb
4099 ! real(kind=8) :: var,ene
4101 open(icsa_in,file=csa_in,status="old",err=100)
4102 read(icsa_in,*) nconf
4103 read(icsa_in,*) jstart,jend
4104 read(icsa_in,*) nstmax
4105 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4106 read(icsa_in,*) nran0,nran1,irr
4107 read(icsa_in,*) nseed
4108 read(icsa_in,*) ntotal,cut1,cut2
4109 read(icsa_in,*) estop
4110 read(icsa_in,*) icmax,irestart
4111 read(icsa_in,*) ntbankm,dele,difcut
4112 read(icsa_in,*) iref,rmscut,pnccut
4113 read(icsa_in,*) ndiff
4120 end subroutine csa_read
4121 !-----------------------------------------------------------------------------
4122 subroutine initial_write
4125 ! implicit real*8 (a-h,o-z)
4126 ! include 'DIMENSIONS'
4127 ! include 'COMMON.CSA'
4128 ! include 'COMMON.BANK'
4129 ! include 'COMMON.IOUNITS'
4131 ! integer :: ntf,ik,iw_pdb
4132 ! real(kind=8) :: var,ene
4134 open(icsa_seed,file=csa_seed,status="unknown")
4135 write(icsa_seed,*) "seed"
4137 #if defined(AIX) || defined(PGI)
4138 open(icsa_history,file=csa_history,status="unknown",&
4141 open(icsa_history,file=csa_history,status="unknown",&
4144 write(icsa_history,*) nconf
4145 write(icsa_history,*) jstart,jend
4146 write(icsa_history,*) nstmax
4147 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4148 write(icsa_history,*) nran0,nran1,irr
4149 write(icsa_history,*) nseed
4150 write(icsa_history,*) ntotal,cut1,cut2
4151 write(icsa_history,*) estop
4152 write(icsa_history,*) icmax,irestart
4153 write(icsa_history,*) ntbankm,dele,difcut
4154 write(icsa_history,*) iref,rmscut,pnccut
4155 write(icsa_history,*) ndiff
4157 write(icsa_history,*)
4160 open(icsa_bank1,file=csa_bank1,status="unknown")
4161 write(icsa_bank1,*) 0
4165 end subroutine initial_write
4166 !-----------------------------------------------------------------------------
4167 subroutine restart_write
4170 ! implicit real*8 (a-h,o-z)
4171 ! include 'DIMENSIONS'
4172 ! include 'COMMON.IOUNITS'
4173 ! include 'COMMON.CSA'
4174 ! include 'COMMON.BANK'
4176 ! integer :: ntf,ik,iw_pdb
4177 ! real(kind=8) :: var,ene
4179 #if defined(AIX) || defined(PGI)
4180 open(icsa_history,file=csa_history,position="append")
4182 open(icsa_history,file=csa_history,access="append")
4184 write(icsa_history,*)
4185 write(icsa_history,*) "This is restart"
4186 write(icsa_history,*)
4187 write(icsa_history,*) nconf
4188 write(icsa_history,*) jstart,jend
4189 write(icsa_history,*) nstmax
4190 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4191 write(icsa_history,*) nran0,nran1,irr
4192 write(icsa_history,*) nseed
4193 write(icsa_history,*) ntotal,cut1,cut2
4194 write(icsa_history,*) estop
4195 write(icsa_history,*) icmax,irestart
4196 write(icsa_history,*) ntbankm,dele,difcut
4197 write(icsa_history,*) iref,rmscut,pnccut
4198 write(icsa_history,*) ndiff
4199 write(icsa_history,*)
4200 write(icsa_history,*) "irestart is: ", irestart
4202 write(icsa_history,*)
4206 end subroutine restart_write
4207 !-----------------------------------------------------------------------------
4209 !-----------------------------------------------------------------------------
4210 subroutine write_pdb(npdb,titelloc,ee)
4212 ! implicit real*8 (a-h,o-z)
4213 ! include 'DIMENSIONS'
4214 ! include 'COMMON.IOUNITS'
4215 character(len=50) :: titelloc1
4216 character*(*) :: titelloc
4217 character(len=3) :: zahl
4218 character(len=5) :: liczba5
4220 integer :: npdb !,ilen
4224 ! real(kind=8) :: var,ene
4228 if (npdb.lt.1000) then
4229 call numstr(npdb,zahl)
4230 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4232 if (npdb.lt.10000) then
4233 write(liczba5,'(i1,i4)') 0,npdb
4235 write(liczba5,'(i5)') npdb
4237 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4239 call pdbout(ee,titelloc1,ipdb)
4242 end subroutine write_pdb
4243 !-----------------------------------------------------------------------------
4245 !-----------------------------------------------------------------------------
4246 subroutine write_thread_summary
4247 ! Thread the sequence through a database of known structures
4248 use control_data, only: refstr
4250 use energy_data, only: n_ene_comp
4252 ! implicit real*8 (a-h,o-z)
4253 ! include 'DIMENSIONS'
4255 use MPI_data !include 'COMMON.INFO'
4258 ! include 'COMMON.CONTROL'
4259 ! include 'COMMON.CHAIN'
4260 ! include 'COMMON.DBASE'
4261 ! include 'COMMON.INTERACT'
4262 ! include 'COMMON.VAR'
4263 ! include 'COMMON.THREAD'
4264 ! include 'COMMON.FFIELD'
4265 ! include 'COMMON.SBRIDGE'
4266 ! include 'COMMON.HEADER'
4267 ! include 'COMMON.NAMES'
4268 ! include 'COMMON.IOUNITS'
4269 ! include 'COMMON.TIME1'
4271 integer,dimension(maxthread) :: ip
4272 real(kind=8),dimension(0:n_ene) :: energia
4274 integer :: i,j,ii,jj,ipj,ik,kk,ist
4275 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4277 write (iout,'(30x,a/)') &
4278 ' *********** Summary threading statistics ************'
4279 write (iout,'(a)') 'Initial energies:'
4280 write (iout,'(a4,2x,a12,14a14,3a8)') &
4281 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4282 'RMSnat','NatCONT','NNCONT','RMS'
4283 ! Energy sort patterns
4288 enet=ener(n_ene-1,ip(i))
4291 if (ener(n_ene-1,ip(j)).lt.enet) then
4293 enet=ener(n_ene-1,ip(j))
4305 ist=nres_base(2,ii)+ipatt(2,i)
4307 energia(i)=ener0(kk,i)
4309 etot=ener0(n_ene_comp+1,i)
4310 rmsnat=ener0(n_ene_comp+2,i)
4311 rms=ener0(n_ene_comp+3,i)
4312 frac=ener0(n_ene_comp+4,i)
4313 frac_nn=ener0(n_ene_comp+5,i)
4316 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4317 i,str_nam(ii),ist+1,&
4318 (energia(print_order(kk)),kk=1,nprint_ene),&
4319 etot,rmsnat,frac,frac_nn,rms
4321 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4322 i,str_nam(ii),ist+1,&
4323 (energia(print_order(kk)),kk=1,nprint_ene),etot
4326 write (iout,'(//a)') 'Final energies:'
4327 write (iout,'(a4,2x,a12,17a14,3a8)') &
4328 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4329 'RMSnat','NatCONT','NNCONT','RMS'
4333 ist=nres_base(2,ii)+ipatt(2,i)
4335 energia(kk)=ener(kk,ik)
4337 etot=ener(n_ene_comp+1,i)
4338 rmsnat=ener(n_ene_comp+2,i)
4339 rms=ener(n_ene_comp+3,i)
4340 frac=ener(n_ene_comp+4,i)
4341 frac_nn=ener(n_ene_comp+5,i)
4342 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4343 i,str_nam(ii),ist+1,&
4344 (energia(print_order(kk)),kk=1,nprint_ene),&
4345 etot,rmsnat,frac,frac_nn,rms
4347 write (iout,'(/a/)') 'IEXAM array:'
4348 write (iout,'(i5)') nexcl
4350 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4352 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4353 'Max. time for threading step ',max_time_for_thread,&
4354 'Average time for threading step: ',ave_time_for_thread
4356 end subroutine write_thread_summary
4357 !-----------------------------------------------------------------------------
4358 subroutine write_stat_thread(ithread,ipattern,ist)
4360 use energy_data, only: n_ene_comp
4362 ! implicit real*8 (a-h,o-z)
4363 ! include "DIMENSIONS"
4364 ! include "COMMON.CONTROL"
4365 ! include "COMMON.IOUNITS"
4366 ! include "COMMON.THREAD"
4367 ! include "COMMON.FFIELD"
4368 ! include "COMMON.DBASE"
4369 ! include "COMMON.NAMES"
4370 real(kind=8),dimension(0:n_ene) :: energia
4372 integer :: ithread,ipattern,ist,i
4373 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4375 #if defined(AIX) || defined(PGI)
4376 open(istat,file=statname,position='append')
4378 open(istat,file=statname,access='append')
4381 energia(i)=ener(i,ithread)
4383 etot=ener(n_ene_comp+1,ithread)
4384 rmsnat=ener(n_ene_comp+2,ithread)
4385 rms=ener(n_ene_comp+3,ithread)
4386 frac=ener(n_ene_comp+4,ithread)
4387 frac_nn=ener(n_ene_comp+5,ithread)
4388 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4389 ithread,str_nam(ipattern),ist+1,&
4390 (energia(print_order(i)),i=1,nprint_ene),&
4391 etot,rmsnat,frac,frac_nn,rms
4394 end subroutine write_stat_thread
4395 !-----------------------------------------------------------------------------
4397 !-----------------------------------------------------------------------------
4398 end module io_config